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Station  I 

PLASMA  THEORY  and  SIMULATION 

A.  FI  LAMENTATION  IN  CHARGED  PARTICLE  BEAMS 

Lee  Buchanan  has  transferred  to  LLL  where  this  work  will  con¬ 
tinue.  No  further  report. 

B.  DRIFT-CYCLOTRON  INSTABILITY  PARTICLE  SIMULATIONS 
Jae  Koo  Lee  (Prof.  C.  K.  Birdsall) 

Work  continuing  on  the  final  report.  No  further  report. 

C.  LOWER  HYBRID  DRIFT  INSTABILITY  SIMULATIONS  USING  ESI  HYBRID  COOE 
Yu-Jiuan  Chen  (Dr.  B.  I.  Cohen  (LLL)  and  Prof.  C.  K.  Birdsall) 

The  linear  properties  of  the  lower-hybrid  drift  instability  were 

studied  using  a  Id  particle-hybrid  simulation,  as  shown  in  the  last  OPR.  The 

model  is  a  slab  with  a  constant  density  gradient;  the  ions  are  unmagnetized 

particles,  shielded  by  the  strongly  magnetized  electrons  through  the  linear 

electron  susceptibility,  y  .  Ions  are  initially  in  a  steady  equilibrium 

e 

state  with  the  ion  diamagnetic  drift  velocity  cancelled  by  the  E*B  drift, 

corresponding  to  electrostatically  confined  ions. 

The  saturation  mechanisms  of  the  instability  are  studied.  In 

2  2 

the  last  QPR,  Fig.  5  shows  the  complex  frequency  versus  v^/v  .  for  w  /w  *  1 , 

E  t i  pe  ce 

mi/me*l600,  Te  =  0,  Ln/Lg  =  0,  Ln/L^.  =  0,  and  the  mode  number  m  =  5  (i  .e.  ,  kx^ 

=  1//2) ,  which  is  approximately  the  most  unstable  mode  for  the  parameters  we 

used.  v_  is  the  ExB  drift  velocity,  v  .  is  defined  as  /T.  /m. ,  and  L  ,  L„ 
t  -  -  1 1  lino 

and  Ly  are  the  scale  lengths  of  the  density,  magnetic  field  and  temperature, 
respectively.  Since  the  growth  rate  y  is  comparable  to  the  wave  frequency, 


2 


the  mode  has  a  wide  band  width  Aco'vy;  and  thus  quasilinear  diffusion  is  one 
of  the  possible  mechanisms  for  the  saturation  of  such  modes.  Allowing  only 
a  single  mode,  we  compared  the  saturation  levels  of  mode  5  with  the  quasi¬ 
linear  theory  as  given  by  R.  C.  Davidson  (Refs.  1  and  2)  under  assumption 
that  the  spectrum  is  strongly  peaked  about  the  fastest  growing  mode,  for 
v  <  v  . .  Figs.  1  and  2  show  that  the  simulation  data  are  in  agreement  with 
the  quasilinear  saturation  level  in  Eq.  48* 


e 

s 


NT. 


1  +  m  /u> 

pe 


2 

ce 


(1) 


where  N  is  the  number  of  ions.  However,  it  is  found  from  the  evidence  of 
oscillations  of  the  wave  energy  at  the  trapping  frequency  uj^.  (Figs.  2a  and 
2b)  and  the  vortex-like  structure  in  the  ion  phase  space  plots  after  satur¬ 
ation  (Figs.  3a  and  3b)  for  v^v^.  *0.57  and  0.85  that  the  end  of  wave 
growth  was  accomplished  by  ion  trapping.  All  of  the  phase  space  pictures 

are  presented  in  the  wave  frame,  x-v  .  *  constant,  where  v  ,  was  calcu- 

ph  ph 

lated  from  linear  theory.  It  is  noted  that  vortices  in  Figs.  3a  and  3b  at 
small  velocities  are  not  due  to  the  multi  beam  instability,  because  those  vor¬ 
tices  would  appear  at  large  velocities. 

Our  explanation  of  why  the  simulations  showed  that  nonlinear  sa¬ 
turation  was  due  to  ion  trapping  but  saturation  levels  agreed  with  the  quasi  1 
ear  theory  is  as  follows.  In  deriving  the  saturation  level  eg  (Eq.  1), 
Davidson  began  with  an  energy  conservation  equation  and  the  only  rea  invo¬ 
cation  of  quasilinear  theory  seems  to  be  the  specification  that  saturation 

occurs  when  the  distribution  function  has  been  "flattened"  around  the  mode 
In  all  of  Davidson's  equations,  there  is  a  factor  of  /2  difference  in  vt.  due 

2 

to  his  definition  of  T.  =m.v  ./ 2. 

i  i  1 1 


Saturation  field  energy  for  the  most  unstable  LHOI  mode 

as  function  of  v,/v...  The  simulation  data  (dots)  are  in 
E  1 1 

agreement  with  the  theory  shown  by  lines  as  discussed  in 


Electric  field  energy  history  plots  for  m./m  -l600.,  u>  /to  =1 

i  e  pe  ce 

Te=«0.,  kX0=1//2. ,  L^/Lg-0 . ,  Ln/LT=0.,  (a)  vE/vt;=0.57;  the  high 
frequencies  are  due  to  the  multibeam  instability  and  low  fre¬ 
quencies  are  due  to  the  LHDI  ;  the  growth  is  that  of  LHDI  as 
discussed  in  previous  QPR's  (using  w  spectrum),  and  (b)  v^/v^.. 
=0.85;  the  multibeaming  effect  is  very  small  (next  page). 


Ion  phase  space  pictures  in  the  LHDt  wave  frames  after  saturatio 
for  m./me=l600. ,  Te=0.,  L^/L^O.,  L^/L^O 

(a)  .=0.57,  and  (b)  v^/vt.*0.85.  The  dots  were  plotted 

one  of  every  5  points. 
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phase  velocity.  Such  flattening  could  in  principle  be  due  to  a  variety  of 
causes  besides  the  usual  quasi  linear  diffusion,  for  example,  trapping.  Fur¬ 
thermore,  the  trapping  frequency  is  given  by 


where 


at_  =  k  ,  /  — ' 
T  m  \  /  mt 


(2) 


[*™s 

k2V 

m 


is  the  electric  potential,  and  k  is  the  wave  number  of  the  most  unstable 

m 

mode.  From  Eq.  1,  we  obtain 


“T 


5/4 


"  “  \45 *1  Wvt. 


(3) 


where 


-HL 


^  A  +  u)2 

pe  ce 


is  the  lower  hybrid  frequency.  From  Ref.  1,  the  corresponding  growth  rate 
and  wave  number  at  maximum  growth  are  expressed  by 


AH 


8  \/2vti 


l  h  ’ 


(4) 


k  =  u).,/v  . 

m  ti 


(5) 


Therefore,  Eq.  3  can  be  rewritten  as 
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When  v£<vtj>  the  trapping  frequency  is  larger  than  the  growth  rate  and  hence 
the  bandwidth  Aw  as  well;  and  the  ion  trapping  will  be  the  saturation  mech¬ 
anism  whenever  the  fastest  growing  mode  is  dominant.  Then  in  Eq.  1  is 
the  saturation  level  due  to  ion  trapping.  Comparison  of  Eq.  6  with  the  ob¬ 
served  trapping  frequencies  (as  in  Fig.  4)  indicates  fairly  good  agreement 
for  cases  in  which  v,_  <  v  .. 

An  estimate  of  the  fluctuation  energy  at  saturation  may  be  made 
through  the  use  of  the  Fowler  thermodynamic  bound  (Ref.  3);  one  possible 
form  of  this  bound  on  the  saturation  level  can  be  obtained  by  assuming  that 
the  system  stabilizes  via  current  relaxation  (meaning  v^  +  0)  as  v^  drives 
the  instability.  However,  this  bound  is  not  applicable  in  our  simulations 
as  our  model  has  assumed  a  constant  density  gradient  and  v^  drift. 

Finally,  multimode  simulations  have  been  done,  i.e.,  all  the 
modes  are  excited  at  the  initial  stage.  Typical  parameters  are  the  same 
as  those  for  the  single  mode  run,  which  were  given  in  the  last  QPR.  The 
fastest  growing  mode  tended  to  reach  the  same  saturation  level  no  matter 
whether  only  a  single  mode  was  kept  or  all  modes  were  included  in  the  simu¬ 
lations.  In  the  frame  of  the  most  unstable  mode,  the  ion  phase  space  plot 

reveals  a  vortex  formation  about  v  =0  after  saturation  for  v_/v  .  =0.57 

x  E  1 1 

as  shown  in  Fig.  5a.  The  length  of  the  plasma  L=2tt,  and  the  number  of 
grids  NG  =  64.  Also  a  dip  appears  at  the  approximate  wave  phase  velocity 
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in  the  ion  distribution  curve  in  Fig.  5b.  As  v^/vt.=0.85,  the  ion  phase 
space  and  the  distribution  curve  are  presented  in  Figs.  6a  and  6b  for  L  = 

2tt  and  NG  =  64.  The  corresponding  ion  phase  space  picture  for  L  =  4m,  and 
NG=128  is  given  in  Fig.  7.  It  shows  that  trapping  still  occurred  for  a 
finer  mode  spacing,  u-cz.  when  Ak  was  reduced  to  half,  i.e.,  Ak^^*  0.07- 
Figs.  5,  6  and  7  show  the  dominant  mode  at  saturation  to  be  X  “u  1//2, 
which  is  the  fastest  growing  mode. 

Our  simulations  exhibited  strong  ion  trapping  at  saturation  in 
all  cases  in  the  low-drift-velocity  regime  with  v^  < v  . .  These  results 
do  not  necessarily  rule  out  quasilinear  diffusion  as  a  possible  cause  for 
saturation  as  the  simulations  had  discrete  modes  (not  a  continuous  wave- 
number  spectrum);  however,  as  each  mode  has  a  large  frequency  bandwidth 
(due  to  large  y) ,  auto-corre lat ion  time  of  the  electric  field  is  signi¬ 
ficantly  reduced.  Trapping  was  also  obtained  by  D.  Winske  and  P.  C.  Liewer 

(Ref.  4)  in  their  2d  particle  simulations  with  v_  larger  than  v  .. 

E  1 1 
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Simulation  (many  mode)  of  LHOI  witn  N =  16384,  <n  /w  =1., 

pe  ce 

m./m  =»  1600.  ,  T  =  0. ,  vc/v  .  *  0.85,  L  /L.  *  0. ,  L  /L_  =  0.  , 

1  e  e  t  1 1  no  n  i 

A^/Ax  *  0. 1414.  Displayed  are  (a)  the  phase  space  in  the 
most  unstable  mode  frame,  and  (b)  the  ion  velocity  distribu 
tion  function  after  saturation.  v^h  is  the  wave  phase 
ve  1  oc  i  ty . 
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0.  NONLINEAR  PERTURBATION  THEORY  OF  THE  LOWER-HYBRID  DRIFT  INSTABILITY 
Yu-Jiuan  Chen  (Dr.  B.  I.  Cohen  and  Prof.  C.  K.  Birdsall) 

The  linear  theory  of  the  lower-hybrid  drift  instability  is  well 
understood  and  has  been  discussed  in  detail  by  Davidson  et  al.  (Ref.  1). 

When  the  amplitude  of  the  wave  is  small  but  finite  after  a  time  of  the 
order  of  the  inverse  growth  rate,  its  further  evolution  will  be  different 
from  the  linear  exponential  growth.  The  nonlinear  dielectric  response  func¬ 
tion  and  an  analysis  of  the  nonl inear  time  evolution  of  a  single  unstable 
mode  are  derived  sel f-cons istently  by  using  perturbation  theory  to  solve 
the  Vlasov  equation  and  the  Poisson  equation.  The  single-mode  approxima¬ 
tion  is  valid  for  the  instability  close  to  the  stability  limit  which  re¬ 
quires  v£<<vti  (the  low  drift  velocity  regime)  for  the  lower  hybrid  drift 
instability.  Modulation  of  the  Langmuir  wave  due  to  the  nonlinearity  has 
been  investigated  extensively  in  the  single-mode  approximation  (Refs.  2 
and  3).  The  nonlinear  evolution  of  drift-cyclotron  and  drift-cone  instab¬ 
ilities  for  plasmas  very  close  to  linear  marginal  stability  have  been  stud¬ 
ied  in  detail  in  both  theory  (Refs.  4,  5,  6  and  7)  and  simulation  (Ref.  8). 

For  simplicity,  we  use  a  one-dimensional  slab  configuration  with 
wave  propagation  in  the  x  direction,  uniform  magnetic  field  in  the  z  direc¬ 
tion,  and  the  density  gradient  in  the  y  direction.  The  ions  are  unmagnet¬ 
ized  as  the  wave  frequency  and  growth  rate  are  much  greater  than  the  ion 
cyclotron  frequency,  and  electrostatically  confined  with  the  ion  diamagne¬ 
tic  drift  cancelled  by  the  E*B  drift,  v^. 

In  Sect.  II,  the  nonlinear  dielectric  function  is  introduced  by 
solving  the  coupled  Vlasov-Poi sson  equations.  Sect.  Ill  is  devoted  to  de- 
reive  the  time  evolution  of  the  lower-hybrid  drift  instability.  The  nonl in- 


-  15  - 


ear  dispersion  relation  obtained  is  used  to  determine  the  field  energy  level 
and  the  frequency  shift  due  to  the  finite  amplitude  of  the  wave.  Finally, 
our  conclusions  are  given  in  Sect.  IV. 

II.  Derivation  of  Nonlinear  Dielectric  Function 


According  to  the  general  theory  of  the  reductive  perturbation 
method,  we  assume  that  the  distribution  functions  F(y,v,t)  and  the  electric 
potential  <j>(x,t)  can  be  expanded  as 


FS(y,y,t) 


Fj(y,y) 


£  enP^(y,v,t)e 
n=1  n 


i  n0 


+  c.c. 


0) 


and 


$u,o  *  23 


n=1 


en<|>n  (x,  t)e‘  n9  +  c.c. 


where 


F^y.y.t)  =  L  eJF^.(y,v,t) 


j=0 


n  =  0,1  ,, 


,(x,t)  =  23  eJ<j>  .(x,t) 

j=0  nj 


n  =  1 ,2,. 


(2) 


(3) 

(M 


9  =  kx  -  tot  (5) 

and  e  =  (£)(  e<{>/T.)«  1.  T  is  the  ion  temperature,  and  k  and  u  are  the  wave 

number  and  frequency  of  a  single  mode.  Quasi  linear  analysis  indicates  the 

2 

current  relaxation  causes  saturation  for  v_  < v  .  (Ref.  9),  where  v  =  T  /M 

E  t i  ti  i 

is  the  ion  thermal  speed.  However,  the  effect  of  the  current  relaxation  is 


•3W  •*  -* 
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small  for  vg<<vtj-  Therefore,  the  density  gradient  and  v^  are  kept  con¬ 
stant  in  our  derivation.  The  distribution  function  Fn ^ (y , v , t )  can  be 
expressed  as 


Fnj  (y c)  "  n0 (y)  fnj  * c)  • 


The  Poisson  equation  of  the  system  is 


-  V2<j>  =  4Trnoe  J dv  (f'  -  fe)  . 


Substituting  Eqs.  (1)  and  (2)  into  Eq.  (6)  yields 


(nk)2$  =  l»7rn  e  f dv  (f*  -  f£)  . 

n  o  J  -  n  n 


Since  the  characteristic  frequency  of  the  lower-hybrid  drift  instability 
is  much  less  than  the  electron  plasma  and  cyclotron  frequencies,  it  is 
assumed  that  electrons  respond  to  the  wave  linearly,  i.e., 


/-e  2- 

dv  f  =  -  x  (nk,nw){nk)  $ 

-  n  Ae  n 


where  xg  is  the  electron  linear  susceptibility.  We  also  assume  a  zero 
plasma  beta  value  (T  0)  to  neglect  electron  resonance  broadening  which 
can  stabilize  the  instability  (Ref.  10).  Using  Eqs.  (3),  (4) ,  and  (9), 
Eq.  (8)  reduces  to 


[1  +  x  (nk.nw) ]  (nk)%  .  =  4irn  e  f  f '  .  dv 

e  nj  o  J  nj 


The  one-dimensional  Vlasov  equation  for  the  ion  distribution  is 


3F  .  3F  e  3d  3F 
3t  3x  M  3x  3v 
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where  M  is  the  ion  mass.  Since  the  density  is  varied  in  the  y  direction 
only,  i.e.,  3Fl/3x  =  0,  £q.  (11)  is  rewritten  as 


3f '  _  e  3d  3f 1 

3t  M  3x  3v 


(12) 


With  the  assumption  of  a  small  perturbation,  the  above  equation  can  be 
integrated  over  the  unperturbed  orbits,  i.e., 


f 


ii  ii_  dt. 

3x  3  v 


03) 


or 


£  ej 

j-i 


f  .  + 

OJ 


i  (f  .  efn9  +  c.c.)] 
n=1  n’J'n  J 


t  « 


j  =  1  -/-*  n=1 

1=0 


(inkp  c"’9  +  c.c.) 

ni 


df 


°d2HZl  +  +  c  \ 

3v  n'  =  1  \  3v  / 


\dt'  .  (14) 


Mote  that  ion  trapping  is  excluded  under  this  assumption.  The  superscript 
i  is  dropped  from  Eq .  (14).  For  j=1,  comparing  the  coefficients  of  exp(ine) 
of  Eq.  (14)  yields 


f 


01 


-  0  . 


(15) 


and 


(16) 


e*10  .  3foo/3v 
M  v-V 


where  V  s  oo/k.  Subs t  i  tut  i  ng  Eq .  (16)  into  Eq .  (10),  then  we  obtain 


[i  +  xe(k»^)j  k2<i>1( 


i*7rn  e 
o 


/3f  / 3 v 
oo 

o  v-V 


k  Xj  (k,u)<j>10  • 


which  yields  the  linear  dielectric  function  D (k , u>)  as 


D(k,uj)  =  1  +  xg(k»w)  +  Xj(k.“)  • 


(The  integral  sub  0  means  integration  over  zero  order  orbits.) 

We  now  proceed  to  obtain  the  second  order  components.  Assuming 
oo=  oj+i6  and  6--0,  coefficients  of  d.c.  terms  of  Eq.  (14)  give 


,  .  .  e  ik  /  *  1C 

02  '  M  "  25  I  '«  3v 


Substituting^  by  using  Eq.  ( 1 6 ) ,  we  get  the  quasi  linear  modification  to  the 
distribution  function 


etj)  ,  /  3  f  /3v\ 

10  3  I  oo  \ 

M  3v  \  (v-V)2/ 


The  component  f^  results  from  the  second  order  terms  of  Eq.  (14)  for  exp(i9) 


e$  3f  /3v 
1 1  #  oo 

M  v-V 


*  — 
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Then,  following  a  similar  method  used  to  derive  Eq.  (17)  or  Eq.  (18),  we 
obtain  from  Eq.  (10) 


D(k,u)d11  =  0 


Coefficients  of  the  second  harmonic  terms  in  Eq.  (14)  give 


i  /e*io 


v-V  3v  \  v  -  V 


3f  /3v\  e<j>-A  3f  /  3v 
oo  \  .  20  oo 


The  first  term  appearing  on  the  right  side  of  Eq.  (23)  is  the  modification 
due  to  the  bare  second  harmonic  oscillation  of  a  single  wave,  and  the  se¬ 
cond  term  represents  its  shielding  effect.  Similarly,  Eqs.  (10)  and  (23) 
yield 


8k  D(2k,2io) 


_  .  .  /3f  /3v\ 

10  f  1  3  /  oo  1 

M  ■'O  v  •  V  3v  \  v-V  / 


For  the  third  order  component  f^>  Eq.  (14)  gives 


*  /.  3V3''  3V3v  ,  ,  3foo/3v  ,  3ftO/3'- 

f12  ‘  M  (+10  - -  ‘  *10  - T  *12  - T  *20  - T 

V  v-V  v-V  v-V  v-V 


Replacing  f^,  fQ^  and  f^  by  Eqs.  (16),  (20)  and  (23),  we  get  likewise 


f12  * 


e*10\  1  32  /3f oo/3v 


M  /  v  -  V  3 v2  \  (v  -  V)2 


,  ,  ,  /  .  ,  / af  /3v 

11  9/1  9/00 


v-V  3v  V  v  -  V  3v  V  v-V 


8k2D(2k 


i  a  /3f  /3v 

1  3/00 


v-V  3v  \  v-V 
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Substituting  Eq.  (26)  into  £q.  (10)  for  < and  using  Eqs .  (2),  (17)  and 
(22)  yields  the  nonlinear  dispersion  relation 


The  first  term  on  the  right  side  is  the  nonlinear  coupling  of  the  potential 
with  the  quasi  linear  perturbat ion.  The  last  two  terms  are  the  nonlinear 
coupling  and  shielding  effect  of  the  second  harmonic  potential  with  the 
fundamental  perturbation.  With  the  definition  of 


W(z) 


-  ,  3f  /3v 

2  f  oo 

Vti  Jo  v-V 


dv 


(28) 


where  z  = V/v^ .  = w/kv^ . ,  the  quasilinear  term  gives 


1  a2/3foo/3v\,  1  d4W(z) 

2  '  ,2v6.  ^ 

1 1 


f  _!_  Ji 

JO  ,2 


-  V  3v  \(v  -  V)‘ 


(29) 


the  bare  second  harmonic  effect  becomes 


_  X  C  _J _ L  I  1  _L  / 3foo/3V\\ 

2  •'O  v-V  3v  \  v  -  V  3v  \  v  -  V  // 


dv 


1  d  W(z) 

I6v^.  dz 
1 1 


(30) 
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and  the  term  associated  with  the  shielded  second  harmonic  oscillation  is 


8k  D (2k, 2m) 


JQ  v  -  V  3v 

. 


3f  /3v 
oo 


32k2X2  v^ . D (2k, 2m)  \  dz2 


d2W(z) 


By  using  Eqs.  (29),  (30)  and  (31),  we  rewrite  Eq.  (27)  as 


O(k,to)ip 


3 


1 _  1  d  W(z)  + 


2k2X2D(2k,2m)  \  dz 


Li2* 


where  ^  =  and  X^  is  the  ion  Debye  length.  The  quasi  linear  term  and 

the  bare  second  harmonic  effect  are  combined  in  the  first  term  of  the  right 
hand  side  in  Eq.  (32). 


III.  The  Evolution  Due  to  the  Nonlinear  Frequency  Shift 


In  this  section,  we  estimate  the  field  energy  at  saturation 
caused  by  a  finite  nonlinear  frequency  shift  by  solving  the  nonlinear  dis¬ 
persion  relation.  As  the  wave  amplitude  is  very  small,  Eq.  (32)  reduces 
to  the  usual  linear  dispersion 


D(k,w)  =  DR(k,w)  +  i D | (k,m)  . 


Let  us  examine  Eq.  (33)  in  the  low  drift  velocity  regime  characterized  by 


|  y/qj  1  «  1  ,  V  and  v^  *  v 


The  dielectric  function  for  the  Maxwellian  ions  is  expressed  as 


hi 


1  U) 


D(k,to)  =  1  +  — - 

hi  k  X  at  -  kv_ 

ce  0  E 


1 


2  ,  2,2  , 
k  ad  lklvti 


The  real  part  of  the  frequency  is  determined  to  zeroth  order  in 


0R(k,ot) 


0) 


1 


0J 


“ce  k  XD  “  -  kvE 


=  0 


The  solution  is 


",  ,  kv_  =  hi 

k2  +  k2  E 


and  the  growth  rate  y  =»  -0  /(30_/ 3tu)  is  given  by 

I  K  O 


y  =  \i  j 


,2,.  2 
k  /k 
tt  m 


d+k2/k2)3  k  \v 

m  m  t  r 


Jlh 


where 


J_  .  1 

2  2  2 

X*  1+u)  /u> 

D  pe  ce 


is  the  wave  number  of  the  most  unstable  mode,  and 


u  . 
Pi 


"ah 


/I  +  hiz  / hi 1 

pe  ce 


(35) 

|  Y/oij  by 

(36) 

(37) 

(38) 

(39) 

m 


is  the  lower  hybrid  frequency  (Ref.  9).  When  the  amplitude  of  the  wave  is 
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small  but  finite,  expanding  Eq.  (32)  around  u>  by  replacing  u  with  u>  + 

o  o 

i(3/5t)  and  using  Eqs.  (35)  and  (36).  we  obtain 

i  [(1  +  ia)  -y]ij>  =  (A+iB)|il»|2<ji  (4l ) 


where 


a  =  (3D./3co)/(3D0/3to)  =  -y/uj 

I  o  R  o  o 


(42) 


A 


Uo~kVE)2  U  ,  1 

8kvc  \3  k2X2DB(2k,2w  ) 
E  v  DR  O 


(43) 


and 


B 


k2'0DR(2k'2“o> 


(44) 


By  using  Eqs.  (36)  and  (39),  we  get 

[k2koV2k-2"o>]  '  ‘  7  • 

k 


Substituting  Eq.  (45)  into  Eqs.  (43)  and  (44),  and  using  Eqs.  (29)  through 
(32),  the  relative  strengths  of  the  nonlinear  contributions  from  the  quasi- 
linear  modification,  the  bare  second  harmonic  oscillation  and  its  shielding 
are  given  as 


Vl.  :  A(2k,2uQ)  b  :  A(2k,2(uo)s 


(46) 


and 

8Q.L.  :  B(2k,2aJo)b  :  8 (2k ,2wq) s 


20  :  -15  :  8 


(47) 


For  the  most  unstable  mode,  k  =  k  ,  the  quasi  linear  modification  on  ti.e  ion 

m 

distribution  function  is  the  dominant  effect.  For  y*r  exp(is)  where  both 
r  and  s  are  real,  Eq.  (42)  becomes 

ar  +  rs  =  -  Ar^  ,  (48) 

and 

r  -  ars  -  yr  =  Br  .  (49) 

Eliminating  r,  we  obtain 


ay  A+aB  2 

2  2  r 
1  +  a  1  +  a 


(50) 


where  the  first  term  is  the  linear  correction  to  the  frequency  in  the 

presence  of  growth,  and  the  second  term  is  the  nonlinear  frequency  shift 

2  2 

which  grows  in  time  with  r  (i.e.,  |e<j>/T.|  ). 

Eliminating  rs ,  we  obtain 


/  Y  B  -  aA 

\1  +a2  1  +a2 


(5D 


Integration  of  Eq.  (51)  yields 
r2  .  r  2  ce2''t/'1^2* 


(52) 


whe  re 
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is  the  field  energy  level  at  saturation.  From  Eqs.  (42)  through  (47),) 
it  is  obvious  that  the  quasi  linear  effect  and  the  shielded  second  har¬ 
monic  oscillation  stabilize  a  single  lower-hybrid  drift  wave,  i.e.,  that 
they  lower  the  saturation  level.  The  nonlinearity  due  to  the  bare  second 

harmonic  oscillation  of  the  wave  raises  the  saturation  level.  If  r  >> 

00 

r(t=0)  =  r  ,  Eq.  (52)  gives 


r2(t) 


2yt 

2  T+^T 
r  e 

o _ 

2 

1+^f  e1+t^ 

r 


and  Eq.  (50)  becomes 


s(t)  = 


A  +  gB 
,  .  2 


.2  1+ct2 

r  e 
o 


■ o  1  +  a2 
1  +—  e 

r 


IV.  Conclusion 

The  nonlinear  dispersion  relation  was  derived.  We  obtained 
the  saturation  field  energy  and  the  nonlinear  frequency  shift  at  satura¬ 
tion  caused  by  the  frequency  shift  only.  With  use  of  Eqs.  (35)  through 
(44),  Eq.  (53)  gives 


16  2k 
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Therefore,  the  nonlinear  frequency  shift  will  stabilize  the  lower-hybrid 

drift  instability  at  small  amplitudes  only  if  k«  k  .  For  the  most  un¬ 
til 

stable  mode,  k  =  k  ,  we  expect  that  other  nonlinear  effects  will  be  the 
m 

dominant  saturation  mechanism  such  as  trapping,  quasi  linear  diffusion,  etc. 
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E.  CONTROL  OF  UNWANTED  BEAMING  INSTABILITIES 
Yu  J i uan  Chen  (Prof.  C.  K.  Birdsall) 


f 

[ 

( 
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No  new  work  this  quarter. 

It  is  planned  to  look  into  similar  instabilities  which  occur 
in  magnetized  plasmas  with  a  Maxwellian  or  other  f(v^)  made  up  of  rings 
in  v^  space  at  t=0.  The  multi -ring  dispersion  relation  will  be  solved 
for  complex  lo  and  real  k.  Instabilities  are  expected  even  with  a  Max¬ 
wellian  distribution.  Simulations  will  be  done  to  find  saturation  lev¬ 
els  and  detailed  ring-ring  interaction  in  order  to  aid  in  finding  means 
of  control  of  these  physical  but  unwanted  instabilities. 


F .  TRANSVERSE  OSCILLATIONS  OF  A  CURRENT  SHEET  IN  A  PLASM*  -  TntCRY 
Alex  Friedman 

We  consider  the  problem  of  on  inf  ini tesfmal ly  thin  current  sheet  of 
ions,  in  o  background  plasma  bounded  by  conducting  -.vails. 
The  configuration  is  illustrated  in  Fig.  I;  the  sheet  current  J°  flows 
in  the  +-y  direction,  the  zero  order  magnetic  field  is  in  the  i-z 
direction  for  x  <  0  and  in  the  -z  direction  for  x  >  0.  and  the  system 


-a  0  x->  a 


z 

I 

Figure  1.  Configuration  of  current  sheet  in  plasma. 


is  bounded  by  metallic  walls  at  x  ■  ta.  We  consider  rigid  displacements 
of  the  sheet  in  the  x  direction  given  by  «  exp(-i(*>t,> . 


The  zero  order  magnetic  field  is  given  by: 


qv05(x). 


Cl) 


where  the  zero  order  current  is  entirely  in  the  y  direction: 
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J°  =  q v05{x)y  (2) 

Tho  first  order  current  due  to  tne  sheet  is  then  given  by: 

±1  -  ~  iu6(x)x  ]  (3) 

(there  is  no  term  orislng  from  i  X  8®  because  is  zero  ot  the 

shee  t ) . 

The  equation  of  motion  of  the  sheet  is: 

=  vonp.os«o*  (4) 

where  0  is  the  cyclotron  frequency  q8/mc.  and  q/m  is  the  sheet's 
charge/mass  ratio:  since  the  sheet  cannot  exert  a  force  upon  itself, 

'*  cancelled  by  (s3/3x)fl°. 

The  analysis  proceeds  in  a  manner  analogous  to  that  of  Ref.  1,  The 
plasma  response  is  given  by: 

Pi'  =  *  B°/c,  (5) 

£'  4  v’  X  B°/c  -  0,  (6) 

which  implies,  since  8°  is  zero  ot  the  sheet,  that  £*  is  zero  ot  the 
sheet,  and  thus  the  equation  of  motion  above  does  not  involve  E',  and 
any  excess  electrons  accompanying  the  sheet  (to  provide  charge 
neutrality)  provide  no  current.  The  field  equations  ore: 

C7  X  8'  =  4K(iphii.) , 


(7) 


-  30  - 


c7XE'  =  itfi1 . 


Combining  eqns.  (6)  ana  (8),  taking  the  cross  product  of  eqn.  (5) 
with  8°,  and  combining  the  resulting  equations,  yields: 

“  U2!'  =  v2  [  tViSk  -  ?X(7X8>) J 

«■  x  C  fu  ~  wb'L  O) 

+■  iuV  X  C(±/nec)xB0]. 

r he  A1 fven  speed  vA  is  defined  os  8°/(4np) 1 /z,  ond  we  hove  used  the 
relation  v#  -  v-j/ne.  The  lost  term  on  the  right  of  equation  (9)  is 
zero,  since  it  is  equal  to  ( i«/nec)d/8x(-jx8°),  ond  jx  is  zero  because 
x'(VXS')  is  zero. 

Using  Eqn.  (3)  for  ond  toking  the  z  component,  yields 

[a2/ax2  ¥  u2/v2  +■  (z/v^Hav^/axja/axjB1 

=  (47T/c)q£¥Q^a2/gx2  (2/vA)(3vA/3x)3/3xj6(x) , 

where  the  subscript  z  of  B]  has  been  omitted.  We  define  a  quantity  b 
equal  to  the  first  order  magnetic  field  everywhere  but  at  the  sheet 


i tsel f  by : 


so  that 


b  *  81  -  (4n/c)qev03(x). 


(12) 

£32/8x2  +■  w2 /v2  +•  (2/vA)(3vA/3x)3/3xJb  »  -(4w/c)q«vrt(w2/v«)6(x) . 


Using  the  relotlon  (2/v. ){3v. /ax)(9b/3x) 


(1/v2)3/9x(v2db/3*)  - 


4 


1 
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d2b/dx.z,  this  simplifies  to: 

b  +•  d/dx[(v2/ij2)dt>/dx.2  ~  -(4n/c)qs v06(x).  ^ 1 


8y  integrating  this  equation  across  the  layer  two  jump  conditions  are 
obtained: 


u?  a"x 


5b  ,  4n 

-  2  =  --  **■,. 


(14) 


Cb]  =  o. 


(15) 


We  specify  a  boundary  condition  corresponding  to  zero  plasma 
densi ty  at  the  wall. 


flx 


=  0 


x=ca 


(16) 


A  particular  solution  of  Eqn.  (13)  is 

bpart  -  -’(2n/c)q*v0(tJ/vA)sin(o|x|/yA). 


(17) 


This  solution  satisfies  the  jump  conditions  at  the  layer  but  not  the 
boundary  conditions  at  the  two  walls,  so  we  add  in  o  constant  x  times 
the  homogeneous  solution  cos(wx/vA),  which  is  symmetric  cbout  x  »  0. 
Then,  using  the  boundary  conditions  (16), 

0b/9x  =  T(2^/c)qev0(ij2/vA)cos(t;x/vA)  -  «(cj/vA)sin(ux/vA),  (18) 


whers 


a  =  -(27T/c)q£v0(u/vA)cot(uaAA). 


(19; 

Keeping  only  the  plasma  contribution  to  A*.  os  noted  after 
Eqn.  (4).  we  use  b  and  not  B1  to  find  the  appropriate  value  of  A1  (*=0) : 

-*J*e  =  =  -(2JT/c2)(qefivo/m)(u/vA)cot(tja/vA) ,  (20) 

or.  more  concisely. 

X  tan  X  -  Y,  (21) 

where 

X  ■  ua/vA,  and  Y  =  2nqev^a/mc2v2  >  0.  (22) 

Note  that  Y  is  real  and  positive.  This  equation  has  an  infinite  n’jnber 
of  discrete,  stable  solutions  for  X.  which  can  be  obtained  in  a 
graphical  manner  by  rewriting  Eqn.  (21)  in  the  form: 

cot  X  =  X/Y  (23) 

and  plotting  both  sides  as  functions  of  X  on  the  sane  axes  (see  Fig.  2): 
it  is  easy  to  prove  that  no  unstable  solutions  exist.  Fur thermore.  for 
large  values  of  Y  the  roots  foil  near  X  -  ±x/2,  t3rt/2,  ....  while  for 
small  Y  the  roots  foil  near  X  -  0  (double  root),  tn,  t2nr,  ... 

It  is  not  difficult  to  extend  the  calculation  to  coses  where  the 
layer  is  not  equidistant  from  two  walls  at  x  ■  -a,  ond  x  -  o2.  Such  on 
equilibrium  requires  the  existance  of  a  uni  form  external  I y  supplied 
magnetic  field  (in  this  sense  it  more  closely  resembles  the 


cyl indricol ly-symmetric  current  layer  than  does  the  centered-sheet 
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l 


X/Y 


Figure  2.  Graph! cal  solution  for  root  of  centered-sheet  pr obletn. 

case);  alternately.  it  may  bo  considered  os  arising  from  the 

introduction  of  tho  conducting  walls  after  the  layer  is  in  place,  thus 

freezing  the  equilibrium  flu*  appropriately.  The  particular  solution 

(17)  is  unchanged,  but  now  the  antisymmetric  homogeneous  solution  .must 

be  included;  we  thus  add  In  or  cos(»x/vA)  +■  0  si n(ux/vA),  and  find 

a  =  -(47t/c)qevDfa)/vA  (24) 

tan(wa]/vA)  +■  tan(u»2/vA) 

where  for  the  normal  mode  frequencies  (but  not  for  the  mode  structure) 
the  value  of  0  is  irrelevant  since  the  field  ot  the  layer  is  the  only 
value  of  Importance.  Defining 

K  *  02/0,,  X  =  ua,/vA,  and  r  =  27tq8v^a,/mc2v2  (25) 

we  find  the  normal  mode  frequencies  to  be  given  by 

X  [  tan  X  tan  KX  ]  /  2  =  Y,  (26) 


i 

\ 

i 

1 

i 

1 

I 

t 

it 


where  Y  and  K  are  real  and  positive,  and  with  no  loss  of  generality  we 


can  choose  0  <  K  <  1 . 
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A  graphical  solution  (for  K  =  Q.B)  appears  In 

1  /  ( tan  X  +  tan  KX) 


X/2Y 


Figure  3.  Graphical  solution  of  uncentered-sheet  problem  for  K  =  0.8. 

Fig.  3.  Due  to  cancellations  of  the  two  different  tangent  terms,  twice 
as  many  roots  are  present  as  in  the  centered-sheet  case.  For  small  Y 
new  roots  at  tn/2.  t3rt/2.  etc.  join  those  at  tn.  ±2n,  ...  For  large  Y 

each  root  splits  into  two,  but  these  remain  close  together  anti  near 
trr/2,  t3rr/2  etc.  when  K  •*  \ .  For  large  Y  (dense  pi  :  ), 

i.e.  when  x  is  near  zero,  the  wave  is  localized  to  one  siti*  of  the  sheet 
or  the  other  depending  upon  which  of  the  "split"  solutions  for  X  Is 
chosen.  The  wave  is  "resonant"  with  either  the  left  cavity  or  the 
right,  and  so  is  of  large  amplitude  (has  a  large  S1)  only  in  that 
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cavity,  while  in  the  other  cavity  the  wove  suffers  destructive 
interference  between  the  particular  and  homogeneous  parts  of  the 
solution.  When  Y  is  small,  a  is  large  and  the  mode  is  not  localized,  os 
expec ted. 

The  author  wishes  to  ack.no.vl  edge  useful  discussion  with  M.  Gerver, 
and  with  0.  Harned,  whose  simulations  of  this  system  are  presented 
elsewhere  in  this  Report. 


1  H,L.  Berk  and  R.N.  Sudan,  "E-!ayer  Precession  in  a  Plasma," 
J,  Plasma  Phys.  6.  413  ( 1971 ; . 
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G.  TRANSVERSE  OSCILLATIONS  OF  A  CURRENT  SHEET  IN  A  PLASMA  —  SIMULATION 
Doug  Haroed  (Alex  Friedman,  Prof.  C.  K.  Birdsal!) 

Our  one-dimensional  quasineutral  hybrid  code  QUAD1  (pre¬ 
vious  QPR)  was  used  to  study  the  oscillations  of  a  thin  sheet  of  ions 
propagating-  through  a  cold  background  plasma.  Conducting  wall  boundary 
conditions  were  used.  The  geometry  of  this  conf i gura t i on  is  shown  in 
Fig.  1  of  the  previous  section.  This  problem  is  analogous  to  the  cylin¬ 
drical  problem  of  the  m  =  0  oscillation  of  a  f  i  e  1  d- revers  i  ng  ion-layer 
in  the  limit  of  infinite  radius. 

Simulations  were  performed  by  placing  a  beam  of  one  cell-width 
in  a  uniform  plasma  and  then  applying  a  small  rigid  perturbation  in  the 
x-direction.  Our  simulations  have  demonstrated  the  stability  of  such  a 
system.  The  measured  layer  oscillation  frequencies  were  found  to  agree 
well  with  analytic  results  derived  in  the  preceding  section.  We  will 

define  x  as  the  half-width  of  the  plasma  slab,  v.  as  the  Alfven  velo- 
w  A 

city,  v  as  the  beam  velocity,  and  M,  as  the  total  mass  of  the  beam.  The 
o  b 

charge- to-mass  ratios  for  the  beam  and  the  background  plasma  have  been 
assumed  to  be  identical.  The  analytic  expression  for  the  oscillation  fre¬ 
quency  of  the  beam, 

X  tan  X  =  Y  (l) 

whe  re 

ux 

v  =  YL  t  o  \ 


Y 


(3) 
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2 

was  solved  for  the  lowest  harmonic  with  a  root-solver  (SOLVER  )  and  is 
plotted  in  Fig.  1.  The  values  obtained  from  our  simulations  are  indi¬ 
cated  on  the  graph. 

The  parameter  Y  on  the  right  side  of  Eq.  1  may  be  expressed 
as  the  ratio  of  the  total  background  plasma  mass  to  the  mass  of  the  beam 


Y 


2x  mn 

W  p 


M 


b 


(M 


There  are  two  limiting  cases. 

When  M  /M,  »1  the  inertia  of  the  beam  is  not  important. 

P  b 

An  initial  perturbation  to  the  right  sends  a  compress i ona 1  Alfven  wave 
to  the  right  wall.  After  the  wave  has  reflected  and  returned  to  the 
center,  the  beam  feels  a  force  to  the  left.  This  motion  sends  a  new 
compressional  Alfven  wave  toward  the  left  wall.  The  beam  continues 
to  oscillate  about  the  center  in  this  manner,  with  a  period  equal  to  the 
time  required  for  an  Alfven  wave  to  traverse  the  slab  twice,  i.e., 

*VA 

<->-—=■  (5) 

X 

w 


An  example  of  this  type  of  motion  is  shown  in  Fig.  2. 

If  the  mass  of  the  beam  is  increased  (or  the  background 
density  decreased)  the  inertia  of  the  beam  slows  its  response  to  the 
wave,  reducing  the  frequency.  For  the  case  Mp/M^ «  1  the  plasma  response 
is  negligible  and  the  motion  of  the  beam  is  that  of  a  simple  harmonic 
oscillator,  with  the  only  force  being  due  to  the  magnetic  pressure  gra¬ 
dient  across  the  beam.  The  equation  of  motion  is 


-7 


(6) 
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Integration  across  the  layer  gives 


-  M 


32C 
b  St2 
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8tt  W  \  (x  -  c)2 
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(7) 


where  BQ  i  s  the  equilibrium  magnetic  field.  Writing  £  =  £eILot  and  per¬ 
forming  a  Taylor  expansion  on  the  right  side  of  Eq.  7»  the  oscillation 
frequency  is  found  to  be 


U) 


B 

o 


(8) 


An  example  of  this  type  of  behavior  is  shown  in  Fig.  3.  It  should  be 
noted  that  the  beam  moves  with  sinusoidal  oscillations,  rather  than  with 
the  triangular  oscillations  which  occur  when  the  effect  is  due  solely  to 
Alfven  wave  reflections  (Fig.  2).  Eqs.  5  and  3  represent  the  high  and 
low  frequency  limits,  respectively,  of  the  lowest  harmonic  of  Eq.  1  (Fig.  1). 
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H.  FIELD  REVERSED  PLASMA  SIMULATIONS,  QUAS I  NEUTRAL ,  in  2d 
Doug  Warned  (Dr.  Alex  Friedman  Prof.  C.  K.  Birdsall) 

We  are  currently  testing  a  two-dimensional  quasineutral  code, 
AQUARIUS  (A  QUasineutral  AlgoRithm  for  Ion  Simulation).  The  purpose  of 
this  code  is  to  study  the  behavior  of  systems  characterized  by  large  ion 
gyroradi  i  and  long  time  scales  (t>fl.\  where  ft.  is  the  ion-cyclotron  fre¬ 
quency).  Examples  of  such  systems  are  field-reversed  mirrors  and  ion  layers. 
Although  these  particular  systems  have  cyclindrical  shape,  cartesian 
coordinates  were  chosen  for  AQUARIUS.  While  cartesian  coordinates  make 
the  application  of  cylindrical  diagnostics  and  boundary  conditions  more 
difficult,  they  avoid  the  problems  of  poor  resolution  at  large  radii  and 
the  singularity  of  the  origin  at  small  radii,  often  associated  with  cylin¬ 
drical  coordinates.  In  our  r-9  code,  there  is  added  limitation  that  a 

Courant  condition  (At<— )  must  be  satisfied  throughout  the  system.  This 

VA 

limitation  governs  the  fineness  of  the  grid  near  the  origin  which  in  turn 
may  force  one  to  pay  a  substantial  penalty  either  in  accuracy  at  large 
radii  or  in  computational  time  (if  small  time  step  is  used  to  allow  a  re¬ 
duction  in  Ax).  Cartesian  coordinates  will  allow  the  treatment  of  a  wide 
variety  of  plasma  configurations,  including  infinite  systems,  for  which 
cylindrical  coordinates  are  not  well  suited  (e.g.,  periodic  boundaries 
cannot  be  applied  in  the  radial  direction).  Cartesian  coordinates  have 
added  advantages  in  that  equations  in  the  field-solver  are  simpler  and 
avoid  difficulties  that  can  arise  in  cylindrical  particle  movers. 

AQUARIUS  is  similar  to  our  one-dimensional  code,  QUAD!  (see 
last  QPR) .  It  is  non-linear,  uses  PIC  techniques  to  advance  the  particles 


and  employs  a  quasineutral  Darwin  field  solver.  As  in  QUAD1 ,  the  electr 
is  advanced  with  the  Darwin  approximation  of  Ampere's  law 


converge. 
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Because  the  successive  application  of  Eqs.  4  and  5  consti¬ 
tutes  two  spatial  derivatives  of  the  electric  field,  some  care  was  re¬ 
quired  to  avoid  the  appearance  of  an  alternating-cell  instability.  Such 
an  instability  will  occur  if  E  and  B  are  determined  using  simple  centered 
differences  to  represent  spatial  derivatives  on  a  single  grid.  To  avoid 
this  problem  we  are  using  interlaced  grids  for  the  electric  and  magnetic 
fields.  Centered  differences  can  then  be  applied  effectively. 

The  code  has  been  tested  on  cold  Alfven  waves.  The  waves  were 
found  to  have  correct  frequencies  and  propagation  speeds  along  the  horizonta 
and  vertical  directions,  as  well  as  at  an  angle  of  45°.  Some  error  has 
been  observed  at  intermediate  angles  due  to  the  squareness  of  the  present 
difference  operators  in  the  code.  We  are  presently  using  four-point  oper¬ 
ators  to  represent  derivatives  for  both  the  cur!  and  gradient  operations. 

2  2 

Such  operators  exhibit  angle  errors  proportional  to  (kAx)  ,  which  become 
severe  at  large  values  of  kAx.  We  would  like  to  reduce  the  angle  errors 
presently  seen  in  AQUARIUS.  While  some  improvement  has  been  obtained  by 
smoothing,  it  may  be  possible  to  further  reduce  the  angle  error  by  the 
implementation  of  a  higher  order  operator. 

Twelve-point  operators  can  be  derived  for  the  gradient  and 
the  curl.  Figure  1  shows  the  grid  points  used  in  the  four  and  twelve- 
point  operators.  Derivatives  in  each  direction  may  be  written  as 


■f. 


+  a(Ui+1,j+3+Ui+l,j-3_Ui-l,j+3'Ui-l,j-3) 


3(Ui+3,j-H+Ui+3,j-!_Ui-3,jt-rUi-3,j-1). 
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where  's  the  effective  wave  vector  produced  by  the  twelve-point 
operator.  It  is  possible  to  compute  the  angle  error  (A^)  of  the  twelve 
point  operator.  The  angle  error  is  defined  as  the  sine  of  the  angle  be¬ 
tween  the  actual  wave  vector,  k,  and  the  finite  difference  approximation 
2  »  fol lowi ng  Ref .  2, 
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The  angle  error  will  be  minimized  when  the  numerator, 
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goes  to  zero.  A  Taylor  expansion  may  be  performed  on  the  sines  and 
cos i nes : 


1 

-  s  i  n 

(X/2) 


cos 


(12a) 


(12b) 


When  these  expansions  are  used  in  Eq.  11,  and  terms  beyond  second  order 
are  dropped,  the  following  expression  is  obtained  to  minimize  the  error: 
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1  -  9a  +  36  =  0  . 


03) 


This  equation  provides  a  means  of  determining  a  and  3  to  give  accuracy  in 

4 

angle  to  (kAx)  for  small  kAx.  The  freedom  to  choose  one  parameter  is  still 

available  and  this  may  be  done  to  effect  some  reduction  in  magnitude  error. 

3  3 

If  the  derivatives  in  the  second-order,  magnitude  error  terms  (3  u/3x  and 
3  2 

3  u/3x3y  )  are  assumed  to  be  comparable,  then  the  optimum  values  of  a  and 
3  are  found  to  be  -.04  and  .0533,  respectively.  Figure  2  provides  a  compar¬ 
ison  between  the  angle  errors  (from  Eqs.  9  and  10)  given  by  this  twelve- 
point  operator  and  the  standard  four-point  operator.  In  Fig.  2  the  mini¬ 
mum  angle  errors,  which  occur  at  an  angle  m/8  relative  to  either  the  x  or 
y  axis,  are  plotted  for  each  operator  (a  =-.04,  B  =  .0533  for  the  twelve- 
point  operator) . 
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|  k  |  A  X 


Maximum  angle  error,  A,  versus  j  k | Ax .  The  dotted  line  is  for  the 
four-point  operator.  The  solid  line  is  for  the  twelve-point  oper 
ator  and  shows  substantial  improvement  at  small  | k | Ax. 
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I.  FIELD  REVERSED  PLASM**  STABILITY; 

LINEARIZED  SIMULATIONS  In  30* 

Alex  Fr 'teaman 

(The  following  is  on  abstract  submitted  as  a  contribution  to  "The 
Physics  of  Mirror  Machines, “  a  hondbooK  under  preparation  at  Lawrence 
Livermore  Laboratory.  The  abstract  is  reproduced  in  its  entirety  hare; 
it  may  appear  in  edited  form  in  the  handbook.) 

The  low  frequency  stability  of  strong  ion  rings  and  ax i symmetric 
field  reversed  mi rror  plasmas  is  being  studied,  using  tune-dependent 
computer  simulation  methods.  The  techniques  employed  are  applicable 
over  a  wide  range  of  parameters,  from  the  large-orbit  field  reversed  ion 
ring  (which  might  be  used  to  confine  cdditional  dense,  less  energetic 
fusion  plasma  as  conceived  of  for  the  Astron  device),  through  FRM 
plasmas  with  Rp/af  ~  5.  reversed-f iel d  theta  pinch  plasmas,  and 
spheromak  plasmas  having  much  smaller  nominal  ion  gyrorodii.  The  effort 
has,  until  recently,  been  concentrated  on  ion  ring  conf igurat i ons  with 
Rp /at-  z  I  and  having  a  dense  cold  background  plasma  component.  Present 
work  i s  on  creating  appropriate  FRM  equilibria  with  Rp/a;  »  2  for 
stabi I  I ty  studi es. 

The  simulation  program  is  a  linearized  three  dimensional  hybrid 
code  called  "RlNGHTBRID"  (FRIEDMAN,  SUDAN,  and  DENAVIT  1*78,1979).  This 
program  model s  an  ion  ring  (represented  by  particles)  in  a  olosmo 
consisting  of  a  cold,  uniform  background  ion  component  and  an 
inertialess  electron  component  of  density  appropriate  for  Iccal 
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quosineutral I ty  (both  inodeled  by  fluid  equations).  Since  th-2  background 
plasma  components  ore  assumed  to  be  pressureless,  the  only  zero  order 
current  is  due  to  the  ring  ion  component.  To  lowest  order  the  r  i  r.g 
particles  (or  the  hot  ions  in  an  FRM  run)  ore  ax i symme t r i c  rings  having 
r,z  coordinates  and  r,8,z  velocities.  First  order  per  turbo c i cns  having 
azimuthal  mode  number  l  a  re  considered,  so  that  each  particle  k  is 
deformed  by  an  inf  ini tes imal  displacement  £ke*p{ i Id) .  Fields  end 
currents  are  represented  by  oxisymmetric  zero  order  ports,  plus  first 
order  par ts  varying  as  exp(iJfl),  all  defined  on  on  Eulerian  mesh  in  the 
r—z  plane.  Since  each  simulation  particle  represents  a  set  of  reel 
particles  lying  on  a  nonaxi symme tr ic  ring,  a  considerable  economy  of 
computation  relative  to  a  nonlinear  3d  code  is  possible.  Coupling 
between  modes  of  different  f,  nonlinear  saturation,  and  other 
large-amplitude  effects  cannot  be  investigated  with  the  linearized 
simulation,  which  serves  largely  as  a  replacement  for  linear  theory, 
since  the  latter  is  difficult  for  such  complicated  configurations. 

A  method  of  generating  quiescent  equilibria  through  the  addition  of 
a  resistive  relaxation  term  (-Ojj3A^/9t)  to  the  zero  order  field  equation 
has  been  developed  (FRIEDMAN  and  5U0AN,  1978).  Because  of  the  chootic 
self-field  betatron  motion  (or  bounce  motion)  of  the  hot  ions,  and  their 
limited  number  in  the  simulations,  only  on  approximate  equilibrium  is 
possible;  the  goal  of  the  method  is  the  minimization  of  fluctuations  in 
the  important  macroscopic  moments  (especially  the  current  density). 

Code  performance  has  been  verified  by  studying  the  normal  modes  of 
the  cold  background  plasma,  which  has  been  represented  both  cs  a  fluid 
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and  by  discrete  particles  centered  in  cells  and  on  cell  edges  (with 
zero-order  velocities  set  to  zero}.  A  consistent  boundary  condition 
arises  from  setting  the  tangential  components  of  the  first-order 
electric  field  to  zero  at  the  walls;  this  necessarily  implies  o  nonzero 
normal  component  of  the  plasma  velocity  at  the  wall  in  this  geometry. 
Further  code  verification  has  included  a  study  of  plosma  return  currents 
across  a  magnetic  field,  in  cylindrical  geometry  and  v/i th  no  center 
conductor  (FRIEDMAN,  SUDAN,  and  DENAVIT,  1979}. 

The  stability  of  inf  I ni tel y-1 ong  current  layer  equilibria  hos  been 
examined  (A.  Friedmon,  In  preporat  ion} .  8ath  stable  and  unstable  kin'* 
and  precessianal  motion  hove  been  observed;  the  unstable  t  —  1  motion  In 
a  radially  decreasing  field  can  be  identified  with  the  MriD  precession 
first  noted  by  BERK  and  SUDAN  (197!},  and  described  in  detail  by 
LOVELACE  (1979}. 

The  stability  of  thin,  weak  rings  (i.e.,  rings  which  ore  not 
encircled  by  field  lines}  hos  been  studied.  Anew  tilting  instability 
of  the  weak  ring  —  plasma  system  has  been  observed  in  the  simulations. 
The  mechonlsm  is  similar  to  that  of  the  kink  instability  of  o  strong 
beam;  however,  unlike  the  strong  ring  case,  the  stability  threshold  is 
dependent  upon  the  background  density.  A  simple  heuristic  analytical 
model  which  contains  the  essence  of  this  mode  has  been  developed 
(A.  Friedman,  in  prepara t ion} . 

Effects  of  ergodic  single-particle  orbits  are  observed  in  many  of 
our  simulations  (FRIEDMAN,  1979};  due  to  the  structural  Instability  of 
particle  t roj ec tor i es,  neighboring  orbits  diverge  exponentially  (and 
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noisily)  with  time.  [n  the  zero  order  motion,  this  effect  I ecds  only  to 
o  loss  of  left-right  mirror  symmetry.  However,  in  the  first  order 
motion  the  associated  "growth"  due  to  s ingl e-par t i cl e  instability  con 
often  be  sufficiently  rapid  os  to  mask  the  collective  modes  which  ore 
the  true  objects  of  study  (the  rote  of  orbit  separation  depends  strongly 
upon  the  details  of  the  equi I ibr lum;  see  FINN,  1979).  Seceuse  of  the 
finite  number  of  simulation  particles  employed,  the  random  phases  of  the 
individual  growing  displacements  cannot  force  macroscopic  moments  to 
be  zero  as  they  would  In  a  true  Vlasov  plosmo  with  an  infinite  number  of 
particles.  Similar  effects  hove  been  observed  in  the  crude  FR.V1 
equilibria  we  have  generated  to  dote;  this  has  been  the  major  impediment 
to  more  rapid  progress  in  this  direction.  This  problem  is  olso  likely 
to  occur  in  future  applications  of  nonlineor  codes  to  problsm.3  of  lineor 
stability  If  certain  "quiet-start"  loadings  ore  used  (i.e.  portlcles  on 
circles  initially). 

The  stability  properties  of  some  field  reversed  ion  ring  equilibria 
hove  been  examined  (A.  Friedman,  in  preparot ion) .  One  moderofely  thick 
ion  ring  with  aspect  ratio  of  order  4:1,  for  which  the  si ngl e-pcr 1 1  cl e 
instability  was  not  excessively  rapid,  has  been  studied  in  detail..  This 
equilibrium  is  stable  to  the  MHO  precession  (t  =  1  radial  mode)  because 
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for  the  possible  existence  of  betatron-resonance  instabilities  {FINN  ond 
5UQAN,  1979)  has  been  observed  for  larger  values  of  8.  The  study  of 
thicker  rings  has  proven  difficult  because  the  singl e-par t ic I e 
instability  masks  the  expected  (slower)  growth. 

The  RlNCA  model,  a  2d3v  (axiaymme t r ic)  nonlinear  ccds,  was  used 
earlier  in  studies  of  ion  ring  formation,  equilibrium,  and  compression 
(FRIEDMAN  et.al.  1977;  MANK0F5KY  et.al.  1978).  The  physics  contained  in 
this  program  Is  similar  to  that  of  the  version  of  the  5UPERLAYER  code 
used  for  mirror  studies,  although  the  algorithms  employed  ore  a i most 
entirely  different.  Companion  runs  using  RInCA  ond  SUPER LAYER  nave 
showed  good  agreement,  thus  verifying  the  performance  of  both  codes. 

Future  applications  of  these  simulations  may  concern  Ian  kinetic 
(e.g.  lass  cone)  instabilities  which  ore  present  in  the  infinite  plosmo 
(8YER5  et.al.,  1978).  However,  the  primary  concern  is  with  gross 
configurational  instabilities  (MHO- 1  ike  kink,  sausage,  etc.)  ond  with 
ion  kinetic  modes  whose  instability  depends  upon  the  detailed  shape  of 
the  field  reversed  equilibrium,  (e.g.  betatron  resonance  or  bounce 
resonance  effects). 


*  This  work  has  been  carried  out  in  col  I aborot ion  with  Profs.  J.  Denovit 
of  Northwestern  University  and  R.  N.  Sudan  of  Cornell  University;  more 
recently,  the  author  has  gained  much  from  interaction  with 
Or.  J.  A.  Byers  of  LLL.  The  author  gratefully  acknowledges  the  support 
of  Prof.  C.  K.  Birdsall  (U.C.  Berkeley). 
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J.  PULSAR  SPARKING 
Or.  W.  M.  Fawley 

A  great  deal  of  progress  was  made  in  simulating  the  initial 
growth,  saturation,  and  final  states  of  electron-positron  pair  creation 
cascades  (= sparking)  in  pulsar  magnetospheres.  In  brief,  it  appears  that 
sufficiently  copious  pair  creation  always  shorts  out  background  electric 
fields  (i.e.,  those  induced  by  the  rotation  of  the  pulsar  magnetic  field) 
and  that  charge  of  both  signs  will  periodically  be  able  to  flow  out  through 
the  light  cylinder.  For  all  initial  conditions  studied  to  date,  the  initial 
burst  of  pair  creation  turns  itself  off  by  shorting  out  electric  fields  all 
the  way  to  the  neutron  star  surface.  The  accelerating  electric  fields  build 
up  again  only  after  the  re  1  a  t  i  v  i  s  t  i  ca  1  I  y  hot  pair  plasma  becomes  sufficiently 
rarified  via  expansion  through  the  light  cylinder  that  the  plasma  suffers 
what  is  essentially  a  dielectric  breakdown.  Due  to  transit  time  effects, 
the  pair  creation  rate  then  oscillates  about  an  equilibrium.  It  is  not  clear 
yet  as  to  whether  these  oscillations  fully  clamp  out  to  zero  or  whether  a 
certain  amplitude  is  maintained  indefinitely.  The  oscillations  have  a  time 
scale  (10  ^  to  10  ^  seconds)  suggestive  of  the  microstructure  observed  in  many 
radio  pulsars.  A  more  detailed  report  will  be  given  in  a  later  Q.PR. 

K.  DIGITAL  FILTERING  IN  TIME 

Dr.  W.  M.  Fawley  (Prof.  C.  K.  Birdsall) 

We  have  made  some  progress  in  developing  time  filtering  algor¬ 
ithms  and  give  a  short  report  here.  A  more  detailed  version  will  appear  in 
a  later  QPR. 
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Following  the  ideas  of  Denavit  et  a!.,  we  have  employed  a  back¬ 
ward-biased  electric  field  in  the  acceleration  subroutine  of  the  1-d  electro¬ 
static  code  ESI.  However,  we  used  a  somewhat  different  biasing  scheme.  Using 


uN+jf-vN~*  q/0 

At  m  ' 


+e)  -N+l  ,  (1-c) 


to 


.  N|+ 1 

requires  a  predictor-corrector  scheme  to  determine  E  which  we  prefer  to 
avoid.  Instead  we  found  that  substituting 


1  ( 


(1+e)  CN+1 /M  ^  (1-e)  _N-1/M' 


(2) 


for  the  right  hand  side  of  Eq .  (1)  for  M>1  reproduced  all  the  useful  damping 
characteristics  of  the  original  equation.  For  sufficiently  large  M  (we  used 

M=5  and  e=0.5  in  most  runs),  it  is  not  necessary  to  use  a  predictor-corrector 

N±1/M  N-i  N  N  N-l 

to  determine  E  inasmuch  as  one  may  advance  V  to  V  using  E  and  E  , 

and  then  fa i r 1 y  accurately  pred i ct  XN±,/,M  and  thus  EN±1/^M.  Computationally, 
this  method  is  somewhat  faster  than  a  pred i ctor-corrector  with  only  one  cor¬ 
rection  step  and  is  much  faster  if  convergence  of  the  predictor-corrector  scheme 
requires  more  than  one  iteration. 

We  have  tested  a  version  of  ESI  employing  Eq.  (2)  on  cold  plasma 
oscillations  and  find  that  for  .jpAt<2,  the  algorithm  is  stable  and  that  the 
damping  rate  increases  quadratical ly  with  oscillation  frequency.  The  algor¬ 
ithm  is  not  generally  stable  for  ii)^At>2,  though  a  "hybrid"  version  using  many 
short  time  steps  followed  by  a  single,  moderately  long  time  step  is  stable  for 
some  parameters  regimes. 


The  damping  algorithm  also  appears  to  inhibit  the  growth  rate 
of  the  multibeam  instability  at  low  spatial  mode  numbers,  but  it  does  not 
suppress  the  instability  at  higher  mode  numbers.  We  plan  to  investigate 
this  problem  further  with  Y.  J.  Chen. 

L.  NONPHYSICAL  BEHAVIOR  OF  HYBRID  OSCILLATIONS  DUE  TO  ALIASING 
Vincent  Thomas  (Prof.  C.  K.  Birdsall) 

While  simulating  hybrid  oscillations  of  a  cold  plasma  excited 
at  large  kAx,  large  variations  in  the  total  energy  were  observed.  These 
variations  were  due  to  an  overall  variation  in  the  kinetic  energy;  there 
was  no  overall  change  in  the  electrostatic  energy  during  the  simulations. 
Moreover,  the  field  energy  was  contained  entirely  in  the  mode  initially 
excited.  These  characteristics  are  shown  in  Figs.  1,  a  through  e.  All  figures 
are  from  simulations  with  the  following  parameters  unless  otherwise  noted: 
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=  14 
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"c  =  3 
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=  E-06 
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The  growth  in  the  kinetic  energy  was  imperceptible  when  the 

lower  modes  were  excited,  but  increased  monoton ica 1 1 y  to  become  many 

times  the  initial  kinetic  energy  when  the  higher  modes  were  excited  initially 

(Fig.  2).  The  frequency  of  the  kinetic  energy  variation  (beat  frequency) 

changed  gradually  from  approximately  0.1  of  the  hybrid  frequency  for  modes 

where  the  growth  has  just  become  noticeable  to  less  than  0.01  of  the  hybrid 

frequency  for  higher  modes.  This  beat  frequency  variation  is  fit  well  by 

the  quantity  w  /2oi  (Fig.  3).  The  plasma  frequency  and  the  hybrid 
P  H 

frequency  are  taken  to  be  the  mode  dependent  quantities  in  this  expression. 


FIG.  I  (a)  Total  electric  field  vs 
time;  mode  II  initially  excited  (b) 
mode  II  electric  field  energy  vs  time, 
(c)  -*■  (e)  kinetic  energy  vs  time.  If 
the  simulation  ran  longer,  the  kinetic 
energy  would  be  seen  to  return  to  its 
original  value. 


•  • 
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The  total  energy  variations  due  to  the  kinetic  energy  variations  were 
much  less  than  the  energy  variations  due  to  the  overemphasis  by  the  grid 
of  the  e lectrostat i c  energy  for  the  first  6  modes  of  a  16  mode  system. 

In  this  sense,  the  kinetic  energy  variation  is  not  important  for  these  modes. 

Phase  space  plots  showed  generation  of  aliasing  (Fig.  4).  Plots 
of  v^  vs.  Vy  displayed  differences  in  phase  that  are  not  to  be  expected 
if  the  motion  is  that  of  anormal  mode  (Fig.  5).  Histories  of  a  single 
particle  velocity  space  motion  showed  alternating  growth  and  decrease  with 
no  drift  (no  E/B  mistake)  (Fig.  6).  The  simulation  run  for  Fig.  6  was  less 
than  the  period  of  an  overall  beat  variation  so  that  the  particle  is  seen 
to  be  still  spiralling  outward.  If  the  simulation  was  run  longer,  then  the 
particle  radius  would  seen  to  return  to  its  original  value.  The  particles 
showed  different  growth  rates,  with  longest  growth  in  positions  where  the 
amplitude  of  the  electric  field  was  largest. 

When  the  time  step  was  changed  from  uj  At=0.15  to  oo  A  t—  1  . 5 »  the 

H  H 

amplitude  and  frequency  of  the  energy  variation  changed  only  very  slightly. 


When  the  mode  number  was  held  constant  and  the  number  of  grid 
points  increased  (meaning  kAx  decreased),  the  variations  decreased,  becoming 
negligible  when  enough  grid  points  were  added.  "Enough"  means  that  kAx<T:/2. 

The  nonphysical  behavior  remained  for  all  excitations  in  x^(0) 
tried  from  0.01  of  a  grid  space  (just  a  little  less  than  the  particle  separa¬ 
tion)  to  10  ^  of  a  grid  space  (which  was  approximately  10  ^  of  the  particle 
separation.  The  frequency  of  the  energy  variation  did  not  change  when  the 
excitation  was  varied  and  the  relative  amplitude  remained  the  same. 


u>  was  varied  from  0.5  to  4  with  the  w  at  1. 
c  P 


Increasing  the 


cyclotron  frequency  yielded  a  larger  amplitude  and  a  longer  period.  Exhaus¬ 


tive  studies  were  not  done,  but  for  all  cases  tried  the  frequency  of  the 
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variation  was  given  by  .  Through  Fourier  analysis  of  the  kinetic  energy, 

a  high  frequency  component  was  also  observed  at  ,  which  can  be  observed 

if  j.r'-u)  is  not  much  smaller  than  w.,  and  so  is  separated  from  the  peak  at 

HC  n 

2<V 

Changing  the  number  of  particles  from  2048  to  256  did  not  change 
the  relative  amplitude  of  the  variation. 

In  an  attempt  to  investigate  possible  beating  between  the  cyclo¬ 
tron  frequency  and  the  hybrid  frequency  (which  would  have  a  large  period  due 
to  the  decrease  in  the  effective  plasma  frequency  at  high  mode  numbers),  a 
simulation  was  done  at  mode  1  with  a  plasma  frequency  corresponding  to  the 
effective  plasma  frequency  of  mode  1 4 .  No  nonphysical  results  were  observed 

which  implies  that  the  poor  aliasing  properties  of  the  higher  modes  is  essen¬ 
tial  for  the  generation  of  the  nonphysical  results. 

One  simulation  was  done  by  perturbing  the  x  velocity  only  but 
the  unphysical  results  still  occurred. 

The  frequencies  of  the  energy  variations  can  be  calculated  by 
making  use  of  the  fact  that  the  electric  field  is  sinusoidal  and  does  not 
change  its  amplitude.  An  essential  feature  of  this  calculation  is  that  it 
depends  upon  the  nons i nusoida 1 i ty  in  space  of  the  force  (or  electric  field) 
at  high  mode  numbers.  A  typical  electric  field  and  its  Fourier  transform 
(treating  the  electric  fields  as  a  continuous  function,  rather  than  a  grid 
quantity)  are  shown  in  Figs.  7  and  8.  For  lower  modes  only  the  principle 
term  would  be  present.  The  analysis  follows. 

The  equations  of  motion  for  hybrid  oscillations,  with  B=zB  ,  are 


(1) 


0 


30 


20 

Position 


Electric  field  at  time  zero.  Note  that  NG  =  32,  mode  14  only  excited 
kinetically  (i.e.,  intial  particle  displacement  and  velocity)  pro¬ 
duces  E  at  mode  14,  plus  other. 


/n|E(J  (arbitrary  units) 
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qE 

x  =  - —  +  uj  y 

m  c" 


Integrating  Eq.  (1)  with  the  special  initial  conditions  of  y  (0)=-uj(_x  (0) 
and  x(0)=0  gives 


y  =  -u)  x 

c 


Putting  Eq.  (3)  into  Eq.  (2)  yields 


••  2  qE 

x  +  a)  x  =  -3— 

c  m 


Fourier  transforming  this  equation  in  x  yields. 


xk  ^c\ 


S(-k)  E (k , t ) 


Here  k  is  to  be  treated  as  a  continuous  variable,  going  from  1  to  “>.  Note  that 

2 

the  r.h.s.  is  not  replaced  by  -to  x^.  From  the  simulations,  we  have  observed 

2  2  2 

that  E(t)  goes  as  cos(«LJt)  where  u)u  =  ou  +  to  (k)  ;  therefore,  we  take 

M  n  c  p 


E  (k ,  t )  =  cos  co^t 


The  most  general  solution  to  Eq .  (5),  using  initial  condition  x^(0)=0,  is 

S(-k)  Ak 

xk(t)  =  C^os  u>ct  -  q  - 2 -  cos  a)Ht  (7) 

mw 

P 


For  those  modes  which  are  not  initially  excited  in  displacement 

or  velocity  (but  are  excited  in  that  an  A^  exists),  X|<(0)=0,  y^COj'O,  we  have 
See  Chap.  8  in  C.  K.  Birdsall  and  A.  8.  Langdon,  Plasma  Physics  via  Computer 


LI 


Simulation. 
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x,  (t)  =- 

k 


Ak  q  S  ( - k ) 


mu 


j^-cos  u^t  +  cos  w^t^ 


(8) 


The  velocity  components  are  vxk=*k  and  Vyk==\=-LlJcxk  w^'c^  produces 

lO  O 

(KE), 


1  /  2  2  , 
'k  2m  Vxk  +  vyk 


-  u 

■\i  2u^_  +  ~2~  ( 1 


cos  2u1Jt)  +  w  ((ui-u  )  cos  (wu+w  )t 
n  C  n  C  n  c 


"  (“u+w  )  cos  (u  -0)  )t]  . 
n  C  n  c 


(9) 


The  high  frequency  at  2wu  is  clearly  seen  and  the  low  (beat)  at  (uu  -  u>  )  x 

H  n  C 

2 

o>p/2uc  is  clear;  the  frequency  (u^+u^)  was  also  observed  through  Fourier 
analysis  of  the  kinetic  energy. 

However,  if  we  are  sufficiently  clever  to  excite  x^  at  t=0  so 
that  C^O,  that  is,  with  proper  normal  mode  amplitude 


x.  (0) 
k 


q  S(-k)  A, 


mo) 


do) 


then  the  solution  is 


xk(t) 


q  s  (-;•■)  E (k) 
2 

mu 

P 


cos 


V 


(11) 


with  no  beating;  only  one  Ak  is  excited. 

The  point  is  that  the  force,  F(x,t),  seen  by  the  particle  (inter 
polated  from  the  grid),  is  never  purely  sinusoidal.  F (k,  t)  =  S  (-k) E  (k)  will 
contain  more  and  more  harmonics  as  the  particle  excitation  is  made  at  larger 
and  larger  kAx-—r.  For  excitation  at  small  k^x«  t  ,  F(k,t)  will  contain  essen¬ 
tially  only  the  excitation  wavenumber  k;  this  will  be 
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q  S(-k)  A 

x  (0)  * - —± 

k  2 

mu 

P 


because 


q  S(-k) 


-  E  (k) 
m 


A  lesson  from  this  exercise  is  that  initial  excitation  at  large 
kAx  is  highly  undesirable,  producing  large  alias  fields  and  subsequent  non¬ 
physical  beat  motion.  In  addition,  as  plasmas  excited  at  small  wavenumber 
will,  due  to  nonlinearities,  excite  modes  at  harmonic  k‘s  and  also  produce 
nonphysical  beats,  it  is  recommended  that  smoothing  be  used,  probably  for 
al 1  modes  kAx  >  r/2. 


r  v 
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Se.c£Lon  II 

CODE  DEVELOPMENT  and  MAINTENANCE 


A. 

ESI 

CODE 

(No  special 

repor  t 

thi  s 

quar  ter) 

B. 

D/ll 

CODE 

(No 

spec i a  1 

repor  t 

thi  s 

quar  ter ) 

C. 

EZOHAR  CODE 

(NO 

special 

repor  t 

this 

quar  ter) 

D.  ES1+-EFL  CODE 

(No  special  report  this  quarter) 

E.  RINCHYBRID  CODE 

Alex  Friedman 

A  new  version  of  the  user's  manual  for  this  code  is  available;  this 

version  contains  a  description  of  the  functions  of  all  code  parameters 

and  variables.  To  obtain  a  copy,  type  (on  the  CDC-7600): 

FT LEM  R05  . CURRENT  MANUAL  /  t  v 

NER-UT  [use]  MANUAL  ULC.  SOX  [r.rn]  VA’ .UAL  /  t  v 

user's  location,  and  [r,nn]  the  box  number. 


where  [use]  is  the 
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F.  RAO  I AL  SIMULATION  CODE  ES1RB  i 

Niels  Otani  (Prof.C.  K.  BirdsaM) 

A  radial  plasma  simulation  code  has  been  developed  for  the  pur¬ 
pose  of  improving  simulation  methods  in  an  r-9  coordinate  system.  Thus  far, 
the  code  has  been  tested  on  simple  plasma  problems.  These  methods  are  in¬ 
tended  to  be  useful  in  the  study  of  a  wide  variety  of  problems  including 
intrinsically  radial  problems  in  both  the  tokamak  and  magnetic  mirror  fu¬ 
sion  machines  and  in  other  cylindrical  devices,  such  as  magnetrons. 

Our  code  is  a  1  J-d  {r,v  ,v  )  code  resembling  in  many  respects 

r  o 

the  cartesian  electrostatic  code  ESI,  developed  by  A.  B.  Langdon  (1970). 

However,  our  code  operates  without  the  grid  used  by  ESI  on  which  charge 
density  is  accumulated  and  fields  are  calculated.  Instead,  particle 
accelerations  are  computed  directly  from  particle-pair  interactions.  The 
main  purpose  for  this  is  to  study  the  effects  of  finite  particle  width 
independently  of  grid  effects. 

Particles  used  in  this  code  are  cylindrical  shells  with  finite 
thickness  with  shape  character! zed  by 

{  —  for  r .  <  r  <  r.  +  w 
w  i  i 

0  otherwise 

where  X  is  the  charge  per  unit  length  of  the  particle  in  the  axial  direction. 

S(r,r.)  is  the  "radial  profile"  of  a  particle  said  to  be  at  coordinate  r. 
and  is  defined  by  S (r , r  )  = 2mrp .  (r) .  The  charge  density  corresponding  to 
this  radial  profile  is  shown  in  Fig.  1.  As  shown  elsewhere  in  this  QPR, 
this  is  the  only  particle  radial  profile  which  will  conserve  both  energy 
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1  Particle  charge  density  as  a  function  of  radius. 


2  Schematic  illustration  of  ES1RB  method  of  advancing  particles 
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3  Summary  of  the  time  centered  leapfrog  method  used  by  ES1RB 

(1)  Jr  -  accel  ,  rotation,  £-accel 

(2)  radial  coordinate  mover 

(3)  velocity  component  mover 
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and  momentum  away  from  the  origin  in  the  limit  At-’-o,  and  therefore  is 
the  natural  choice  for  the  particle  shape. 

The  equations  of  motion  to  be  simulated  are 
2 

Va  1  f 

<!  =  —  +  -  /  S (r)  E  (r)  dr  +  a  v  (la) 

r  r  m  J  r  c  B 


V9Vr 


U)  v 

c  r 


Ob) 


where,  by  Gauss'  Law, 


E(r) 


2 

r 


E /^(r' 

j  Jo 


dr' 


m  is  the  mass  per  unit  length  in  the  axial  direction,  and  u>c  is  the  cyclotron 
frequency,  AfWmc.  In  other  words,  we  consider  only  forces  derived  from 
a  uniform  external  axial  magnetic  field,  and  from  a  self-consistent  elec¬ 
tric  field,  as  in  ESI . 

The  ACCEL  routine  of  our  code  employs  Langdon's  ha  1 f-accel era- 
t ion-rotation-ha  1 f-acce lerat ion  scheme  modified  for  cylindrical  coordinates: 


Vr1  = 

<\>old  + 

1  41 

m  2  • 

fsM  tr 

(r)  dr 

II 

CM 

L_ 

> 
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(half-acceleration)  (2a) 


(rotation)  (2b) 


(half-acceleration)  (2c) 


A  time-centered  leapfrog  method  of  advancing  particles  is  used: 


(v)  u'vft),  (v)  =  v(t+At),  and  r  =  r(t+At/2).  Provisions  also  exist  in 

-old-  -  new 

this  portion  of  the  code  for  adding  velocity-proportional  damping  and  uni¬ 
form  background  charge. 

The  MOVER  advances  both  the  coordinate  and  the  velocities  to 
the  new  position  according  to  the  formulae: 


r  +  ~  = 

new  2 

y(roid+  %  + 

vrAt)2  +  (vgAt)2 

(3a) 

(vr) new  = 

(vr)old  COS  9 

+  ^old  sin  9 

(3b) 

(v0 } new  = 

"  (vr) o 1 d  Sin 

9  +  (Vold  COS  6 

(3c) 

where 


cos  9 


v0At 


r  +  w/2 
new 


s  i  n  9  = 


r  +  w/2  +  v.  At 
old  9 


r  +  w/2 
new 


(3d) 


(3e) 


These  relations  are  derived  from  the  schemati c  i  1 1 ustrated  in  Fig.  2. 
We  wish  to  have  the  center  of  mass  of  each  angular  section  of  the  cylindrical 
particle  moving  according  to  Eq.  (1).  As  is  evident  from  Fig.  2,  the  center 
of  mass  coordinate,  r+w/2,  changes  according  to  (3a),  and  the  components  of 
velocity  vector  v  are  rotated  by  angle  9,  owing  to  the  rotation  of  the  coordi¬ 


nate  system  at  the  new  position  relative  to  the  old.  Note  that  the  rotation 

does  not  advance  v  in  time.  It  is  however  equivalent  to  the  -v,v  /r  term  in 

6  r 
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Eq.  (1b),  and  the  centrifugal  force  term  in  Eq .  (la).  The  steps  involved  in 
advancing  the  particles  are  summarized  in  Fig.  3- 

Initialization  of  the  particles  is  achieved  by  means  of  an  input 
file  containing  r,  vf,  and  v^  values  of  all  the  particles.  Pre-processor 
programs  have  been  developed  to  generate  these  input  files.  This  method  of 
loading  particles  allows  a  great  deal  of  flexibility  with  a  minimum  amount 
of  effort.  For  instance,  a  simple  method  for  establishing  a  one-dimensional 
radial  equilibrium  is  to  load  the  particles  arbitrarily  and  use  ES1RB  itself 
with  damping  to  damp  the  particles  to  equilibrium.  ES1RB  then  produces  a 
record  of  the  equilibrium  positions  and  this  is  used  as  the  input  file  for 
the  pre-processor,  wh i ch  can  mod i fy  the  equilibrium  in  whatever  way  desired. 

The  modified  file  is  then  accepted  as  an  input  file  by  ES1RB  for  the  main  run. 

To  start  the  leapfrog  scheme  a  routine  analogous  to  the  ESI  SETV 
is  called  once  immediately  after  the  particles  loaded.  This  routine  performs 
an  acceleration  followed  by  a  rotation  backwards  in  time  by  an  amount  At/2, 
and  thus  provides  the  proper  offset  between  position  and  velocities. 

Some  Resul ts 

The  equilibria  for  all  simulations  run  so  far  have  been  established 
by  the  damping  method  just  described.  This  process  is  illustrated  in  Fig.  4. 

It  was  found  that  the  desired  equilibrium  was  most  efficiently  achieved  using 
a  two-stage  damping  scheme.  The  first  stage  used  a  decay  constant  of  m  (the 

plasma  frequency)  while  in  the  second  stage,  damping  is  reduced  to  O.lu  . 

P 

Radial  plasma  osci a  I lati ons  were  simulated  quite  nicely  by  esta¬ 
blishing  equilibrium  in  the  presence  of  a  uniform  background  charge  distribu¬ 
tion.  This  equilibrium  was  perturbed  slightly  and  used  as  the  initial  configur- 


[•It] 


[•] 


Energy  Radial  coordinate 

Equilibrium  for  the  main  run  is  established  using  this  preliminary 

damping  stage.  Only  one  stage  was  needed  for  thie  particular  run. 

Number  of  particles  = 50,  2ir/u  =4.42,  damping  constant  =  1  .0w  , 

P  P 

width  of  parti  cl es  = 0. 05 . 
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ation  for  the  main  run.  The  result  is  shown  in  Fig.  5-  The  Fourier  spec¬ 
trum  of  the  kinetic  energy  peaks  quite  sharply  at  2ui  ,  as  expected.  We  also 

P 

find  that  the  use  of  wider  particles  and  larger  perturbations  produces  a 
more  complicated  kinetic  energy  spectrum,  for  reasons  not  yet  completely 
understood . 


Using  a  similar  method,  hybrid  oscillations  were  also  simulated 

(Fig.  6).  The  kinetic  energy  spectrum  shows  components  at  as  well  as 

2oj  .  This  was  found  to  be  a  result  of  the  initialization  of  all  particles 
P 

with  v,  =0.  This  leads  to  a  superposition  of  a  E*B  drift  on  the  hybrid  os- 
cillation  motion  of  the  particles.  A  calculation  taking  this  E  x  B  drift 
into  account  quantitatively  agrees  with  the  relative  amplitudes  of  the  j 

P 

and  Fourier  peaks  of  4  to  1  found  computationally. 

We  have  also  found  that  the  code  accurately  reproduces  the  equil¬ 
ibrium  configurations  of  a  completely  non-neutral  plasma  column  in  the  pre¬ 
sence  of  an  axial  magnetic  field  as  described  in  Davidson's  Theory  of  Non- 
Neutral  Plasmas,  Ch.  1. 

At  the  present  time  a  study  is  being  made  of  the  possibility  of 
the  presence  of  a  two-stream  instability  in  the  magnetic  insulation  sheath 
surrounding  a  cylindrical  cathode.1 

Of  considerable  concern  in  this  project  is  the  fact  that  an  ade¬ 
quate  method  of  dealing  with  the  passage  of  particles  through  the  origin  has 
not  yet  been  developed.  As  suggested  in  another  article  in  this  QPR,  accurate 
(i.e.,  energy  conserving)  methods  of  simulating  particle  motion  through  the 


origin  may  not  be  easy  to  come  by. 


SUMMARY 
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Work  on  the  problems  involved  in  radial  plasma  simulation  is  con¬ 
tinuing.  So  far  we  have  limited  projects  to  the  simulation  of  known  results 
and  have  met  with  moderate  success.  Some  difficulties  involving  the  passage 
of  particles  through  the  origin  have  been  encountered.  In  the  future,  the 

code  will  be  expanded  to  either  a  gridded  2-d,  r-S  code,  or  to  a  1?-d  gridded 

2 

code  employing  a  Fourier  transform  method  in  3  to  simulate  a  2-d  code. 
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G.  CONDITION  ON  PARTICLE  SHAPES  IN  ID  RADIAL  CODES 
Niels  Otani  (Prof.  C.  K.  Birdsall) 


A  simple  analysis  shows  that,  away  from  the  origin,  only  parti¬ 
cles  with  S (r , r. )=S (r-r . )  can  simultaneously  conserve  energy  and  momentum 
in  the  limit  At-H).  (Here  S(r,r.)  is  defined  by  S(r,r^)=2irp.  (r)  for  a  parti¬ 
cle  located  at  r..)  If  particles  are  forced  to  pass  through  the  origin, 
energy  conservation  is  impossible. 

For  energy  to  be  conserved,  we  require 


dK  dU 
dt  dt 


0 


(1) 


where  K  is  the  total  kinetic  energy: 

K  =  ]£  }mv? 

i 


(2) 


in  which  v.  is  the  radial  velocity  of  the  i-th  particle,  and  U  is  the 
total  electrostatic  potential  energy: 


U 


27rrdr 


E*(r) 


8tt 


(3) 


It  is  understood  that  we  are  working  in  cylindrical  coordinates  and  that  all 
extrinsic  quantities  are  per  unit  length  in  the  axial  direction.  We  as¬ 
sume  that  the  electric  field  is  due  to  the  presence  of  all  the  charged 
particles  and  to  a  fixed  background  charge  distribution.  Thus,  from 
Gauss'  Law, 


yn 


f  Vr)  * 


r.)  dr' 
J 


w 
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where  Q.(r)  is  the  background  charge  inside  radius  r.  We  find  then  that 


dK  dv. 

—  =  £mv  — 1  =  £v  f  dr  S(r,r  )  E  (r) 
At  ;  1  At  ;  1  •'n  \  r 


•  f0ir  (E'i  s(r'ri>)  7(Vr)  *  -r  iUr') 


while 


=  1  f  r  dr  2E  (r)  !!r 

dt  kJn  r  dt 


■  l  dr^(x  ^  7^b(r)  */?s<r'’rj)  dr') 


dri  /-r  3S 


*  /  dr(Z  -7T  /  —  (r‘,r  )  dr ’ j  — ((Q.  (r)  +  /"  £s(r',r.)  dr 

•/0  \i  '0  3r.  '  r\  J0  j  J 


'0  3r . 


For  all  practical  purposes,  this  necessitates 


dr.  -r  . 

— 1  /  -  (r',r.)  dr'  =  -  v.  S(r,r  ) 

dt  3r.  ' 


using  Eqs.  (5)  and  (6)  in  Eq.  (l).  In  the  most  general  case,  we  need  not 


have  dr./dt=v!,  indeed,  a  mover  using  an  algorithm  of  the  form 


^  ( ( r  i )  _.  )  =  f((r.)  ,  .)  +  v.At 

i  new  i  old  i 


is  conceivable.  In  the  limit  At-<-0,  we  then  have  v.  «- - .  Allowing 

1  dr.  dt 
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for  this  possibility,  we  see  Eq.  (7)  is  equivalent  to 


3S  ,  , 

”(r’ri) 

I 


S ( r=0 ,  r . ) 


df  3S,  , 

dTT  a7(r’ri) 

i 


If  f(r.)=r.,  then  (8a)  is  equivalent  to 


S(r, r. )  =  S (r-r . )  . 


In  other  words,  if  energy  is  to  be  conserved,  the  radial  profile  S  must  be 
rigid  and  cannot  change  shape  as  the  particle  moves  in  and  out.  Clearly, 
such  a  restriction  is  incompatible  with  Eq.  (8b)  when  particles  pass 
through  the  origin;  therefore  energy  conservation  is  impossible  in  any 
1-d  radial  system  in  which  particles  pass  through  the  origin. 

A  similar  calculation  can  be  done  in  the  spherical  case;  here 
we  find  the  same  theorem  true  when  the  radial  profile  S  is  defined  by 


S  (r,  r . )  =  4it  r  p.(r) 
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H.  RAO  I AL  CODE  NOTES  (R,  R9 ,  RZ,  R9Z) 

1.  INTRODUCTION  (C.  K.  Birdsall) 

Simulations  using  radial  coordinates  or  grids  pose  a  number  of 
problems  for  which  some  solutions  will  be  presented.  The  problems  and  solu¬ 
tions  to  be  presented  are  not  claimed  to  be  new  or  unique  or  exhaustive. 

The  object  is  to  gather  together  in  this  and  succeeding  QPR’s  some  of  the 
problems,  solutions,  and  experiences  with  radial  simulations. 

Radial  coordinates  may  be  preferable  to  rectangular  in  Id  simula¬ 
tions  using  cylindrical  and  spherical  R  only,  in  2d  (cylindrical  with  R9  or 
RZ,  and  spherical  with  R0 ,  R$,  9<t> )  and  in  3d  (cylindrical  with  R9Z,  spherical 
with  RB^) . 

An  object  of  using  radial  coordinates  or  grids  is  to  emphasize 
radial  behavior,  usually  implying  that  the  physical  and  simulation  models 
have  a def i n i te  origin  (R=0) .  There  may  be  strong  radial  forces  and  motion  or 
radial  (circular  or  spherical)  boundary  conditions  such  that  circular  or 
spherical  harmonics  are  more  easily  identified.  Radial  coordinate  use  is 
also  promoted  by  the  concern  that  using  a  rectangular  mesh  with  rectangular 
boundaries  will  introduce  unwanted  moments  in  the  grid  quantities.  This 
argument  may  hold  for  XY  contrasted  with  R0 ,  but  not  for  XY  contrasted  with 
RZ.  However,  some  conventional  wisdom  argues  against  abandoning  rectangular 
meshes,  pointing  out  that  with  XY  meshes,  problems  with  the  origin  are  obvi¬ 
ated,  the  equations  of  motion  are  more  easily  handled,  the  radial  boundary 
conditions  can  be  handled  to  good  approximation  and  the  circular  harmonics 
diagnostics  are  obtainable  to  good  approximation.  Preferences,  we  think, 
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will  be  made  easier  and  more  rational  by  displaying  the  details  of  using  ra¬ 
dial  coordinates  and  meshes. 

The  QPR  radial  code  reporting  comes  from  several  sources,  identi¬ 
fied  section  by  section. 

* 

2.  ONE  DIMENSIONAL  MODELS  (C.  K.  Sirdsall) 


One  dimensional  radial  coordinates  or  meshes  are  clearly  of  gen¬ 
eral  interest  and  eminently  useful.  An  example  is  the  set  of  electrostatic 
plasma  diode  problems  solved  by  Barnes.^  Let  all  of  the  charges  be  cylin¬ 
drical  shells  (spherical  shells  can  be  done  by  inference)  which  move  radially 
due  to  radial  electric  forces.  The  electric  field  within  the  sth  shell  at 
rs  (see  Fig.  1)  may  be  obtained  directly  from  Gauss'  Law, 


The  electric  force  on  the  sth  shell  is  dependent  only  on  the  net  charge  within 

that  shell.  The  radial  force  on  the  s*"^  shell  of  uniform  charge  density  p, 

2  2 

line  density  p^=Trp(b  -a  )  due  to  all  other  charges  is 


F 

r 


dl/  . 


(2) 


Specifically,  the  radial  force  on  Shell  I  due  to  Shell  II  is  (for 
uniform  density  shells) 

Frl  -[°I  ErII  dV  "  /  PI  2wr  dr  dz 


Christopher  W.  Barnes,  "The  Computer  Simulation  of  a  Spherically  Symmetric  Plasma", 
SUIPR  Report  No.  3^.  Inst,  for  Plasma  Research,  Stanford  Univ.,  March  1970. 
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p£i  pm 

2tt£ 


using 


Pt  (b  -  a) 


p j(b  -  a) (b  +  a)  u 

-v  (b  +  a)  _ 

2  2  T 


where 


(3) 


_  a  + b 


That  is,  when  Shell  II  is  wholly  inside  of  Shell  I,  then  the  force  due  to  Shell 
II  acts  as  if  applied  at  the  geometric  middle  of  Shell  I,  independent  of  shell 
thi ckness. 

The  self  force  on  Shell  I,  due  to  Shell  I  (after  some  algebra),  is 


w  5  f 


PI  Erl  dV 


'  2.1  /  a  +  t/3 


dz 


2irer  \  2r 


dz  for  t  «  a  ^  b  'v  r  (M 
s 


This  force  is  always  radially  outward;  this  effect  is  not  found  in  planar  slab 
mode  1  s . 

Let  the  shell  be  placed  in  a  uniform  background  density  of  value 
p^,  opposite  in  sign,  so  as  to  oppose  the  self-force,  with  inward  force  on  a 
thin  shel 1  of 


F.  =  P. 

i  n  l 


dz  . 


(5) 
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The  self  and  background  forces  compete,  as  shown  in  Fig.  2,  going  to  zero  at 
vrequi  librium’  obtained  from 


dz  =  —  ~ 


1  /  pl 


tr  rz  p.  d2  = 


2pl  dz 


L 


(background  charge  enclosed)  =  (half  of  the  shell  charge) 


The  net  force  is 


1  0 o  p0pkr 

—  — _i_+  g-  b  s 

2  2 ner  2s 

s 


2s  pl  pb  (  —  +  r 


P5.pb  (rs  +  re)(rs‘re) 


which  reads 


Several  methods  for  assigning  particles  to  a  radial  grid  will  be 
given  here.  In  all  methods,  the  particles  have  coordinates  r.  and  charges 
q.  which  are  to  be  assigned  to  nearby  grid  points  r. .  The  particle  shapes 
will  fall  out,  with  some  differences  from  rectangular  particles  and  grids; 
e.g.,  the  shapes  (ordensity  profiles)  may  not  be  symmteric  about  a  center, 
or,  the  average  charge  density  of  a  particle  will  decrease  as  the  charge  (of 
fixed  radial  thickness)  moves  radially  outward,  or,  the  assignment  at  the 
origin  may  differ  from  that  elsewhere,  etc. 

METHOD  A  (from  B.  I.  Cohen,  historical  origin  obscure) 

The  charge  q.  is  located  at  r.  between  r.  and  r..,,  as 

3  i  i  j  j+1 
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charge  density,  which  is 


p.  H  q./( volume  of  charge) 

which,  in  cylindrical  coordinates,  for  a  cylindrical  shell,  is 


q. 

(2irr.  5r<$z) 

i 


where  5r,5z  are  the  nominal  fixed  dimensions  of  the  charge.  The  density  to  be 
associated  with  r.  is  linearly  weighted  from  r. ,  as 


where  Ar  =  r . , . 

j  +  1 


•r.  and  the  density  associated  with  r.At 
J  J  +  1 


However,  charge  is  to  be  assigned,  as  follows: 
q(r  J  =  p(cj)  2irr.  5r  5z  , 

q(rj+l)  =  p(rj  +  l)2,T  'j  +  1 


6r  5z 
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which  is 


The  check  on  charge  conservation 
q(r.)  +q(rj  +  1)  =  q. 

is  observed.  Applying  these  weights  produces  q(tj)  as  a  function  of  r.  as 
shown  in  Tig.  1. 

A  similar  approach  works  in  spherical  coordinates  with  spherical 
she  11s,  where 


(4irr?<5  r) 


with  assignment  of  density  done  quadrati ca 1 1 y ,  as 


leading  to 


F/G  20/9 
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FIG.  1 


Charge  assigned  to  r.  as  a  function  of  charge  position  r. .  For 
r.  »Ar,  the  familiar  planar  linearly  weighted  triangular  shape  is 
seen.  For  r.  a  few  Ar,  the  shape  is  distorted,  as  shown. 
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which  also  conserves  charge,  q(r^)  +q(rj  +  1)  =q 

I.  RJET  DEVELOPMENT 

(No  special  report  this  quarter) 
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SOFTWARE  DEVELOPMENTS 

H.  Stephen  Au-Yeung 


(1)  FREX 


FREX  allows  the  user  to  extract  selected  frames  from  one  or  more 
FRB0  graphic  files  on  the  CRAY  into  one  FRH0  output  file. 

This  program  can  be  obtained  by  typing! 
rfilem  read  Iddd  .cray  frextESClena  /  t  v 


This  document  corresponds  to  the  FREX  version  of  October  10. 
ly/y.  Later  versions  of  FREX  will  be  stored  in  FILEIi  directory 
.cray  of  user  number  Iddd.  The  user  should  periodially  check  the 
date  of  this  file  tfilem  how  Iddd  .cray  froxl  to  see  if  the 
program  has  been  updated.  The  file  FREX  is  a  LIB  file!  it 
contains  the  latest  sources  as  well  as  the  latest  documentation. 
To  get  the  documentation,  type! 

lib  frextLFlx  f rex/doc  ILF lend  /  tv 

net out  fuse]  frex^doc  Lswp.3  box  nnn  frex  /  tv 


The  commands  in  FREX  are  similar  to  those  in  DDEX.  a  DDU0  file 
frame  extracting  program  that  resides  on  the  C&C-/b00  tsee  LIBRIS 
*5/1.  The  following  is  a  summary  of  all  FREX  commands! 

L  integer/  - 

The  frame  (.represented  by  the  number  L  integer/1  to  be 
extracted.  IFREX  allows  only  the  second  one  of  the  two 
header  frames  to  be  extracted  by  the  user!  it  is  referred 
to  as  frame  number  0.  The  first  graphic  frame  is  referred 
to  as  frame  number  1.  and  so  on.l 

<.  integer/  thru  <.  integer/  - 

"thru"  can  be  replaced  by  "l"  or  "to".  It  causes  all  frames 
between  the  two  integers,  inclusive,  to  be  extracted. 

infile  Lfile  name/  - 

Change  the  input  file.  Hi  1  selected  frames  of  the  previously 
opened  file  are  extracted  before  opening  the  new  one. 


end 

Terminate,  processing  all  selected  frames, 
quit  - 

Like  "end",  but  when  entered  after  '.he  "infile"  command  causes 
previously  entered  frames  to  be  extracted  without  changing 
the  input  file  or  causing  termination!  more  frame  numbers 
can  then  be  entered. 

box  Lnnn/  Lid/  - 

The  id  string  consists  of  up  to  dA  alphanumeric  characters! 
embedded  spaces  may  appear  incorrectly  ond  are  not 
recommended!  use  the  poriod  or  slash  instead.  The  output 
file  name  is  constructed  from  the  id  lino!  it  is  of 


the  form  fkK  first  six  characters  of  id/.  If  iti  consists 
of  fewer  than  six  characters,  the  filename  is  right-filled 
with  asterisks  to  ensure  an  o ig th-charac lor  file  name. 

If  the  file  to  be  generated  exists,  the  user  is  given 
the  option  of  overwriting  it. 

s ize  .  t  intger >  - 

To  specify  the  size  fin  words!  of  the  output  file  fdefaults 
to  zlUyU UUb!  .  this  command  lias  to  be  issued  before  the 
box  and  id  is  entered!  otherwise,  the  default  size  will 
bo  used . 

family  \ the  first  file  name  of  a  file  family/  - 

This  is  the  same  as  the  “ir.t'ile"  command  except  that  when  the 
end  of  file  is  encountered,  lvRE;<  will  open  the  next  file 
within  the  family.  All  file  names  must  bo  exactly  b 
characters  long.  The  last  name  in  the  family  should  end 
with  an  "x“  fo.g.  f  lUSrpilx  or  fytestdx,  etc.!.  The  last 
two  characters  of  other  files  must  be  numeric  between 
"UU"  and  "yy“. 

Note:  The  BASEL  IB  routine  ZSEOHSI-'’.'  is  used. 

cancel  integer >  - 

Cancel  the  frame  number  < integer/  that  has  previously 

entered  for  extraction  '.frame  must  be  cancelled  1-by-l!. 

nochar  - 

Eliminate  all  alphanumeric  characters  from  all  frames.  Useful 
for  suppressing  crowded  labels  on  axes  when  preparing 
figures  for  publication. 

char  - 

Reverse  the  effect  of  "nochar"  li.e.  return  to  default  mode!. 

offset  <  integer/  - 

Set  offset  for  calculation  of  frame  numbers  Idefaults  to 
13  and  usually  not  required!.  For  example,  if  the  number 
i  is  entered,  all  frame  numbers  entered  after  this 
command  will  have  i  added  to  thorn  before  processing. 


Restrict  ions: 

11J  The  routine  will  fail  if  the  output  exceeds  the  size 
specified  by  the  "size"  command. 


Example  1! 

This  example  shows  how  to  extract  frames  from  different  FRbB 

f  i  1  v.*'.. 


user:  frex  size  400000b  box  b44  solv  /  1  .4 

rout ine/user :  /infile  fl04.a.l 
rout ine/user :  >1 

rout ine/user :  /infile  f 104. a. 4 
routine  :  Processing  file:  f 104.0.1 

routine/user:  >1  4  4 

rout ine/user :  /infile  fl04.b.l 
routine  :  Processing  file:  f 104.0.4 

rout ine/user :  >114 

rout  ine/user :  >end 

routine  :  Processing  file:  f!04.b.l 

routine  :  all  done 

The  output  file  in  this  case  will  be  "f0solv**“. 


Example  4: 

One  thing  not  mentioned  above  is  that  the  first  input  file  appearing 
on  the  execution  line  right  after  "frex"  needs  no  "infile"  command. 

user:  frex  f 105rp0x  box  b44  easel  /  1  .1 
routine/user:  >b  t  40 

routine/user:  /cancel  IS  cancel  41  end 

routine  :  Processing  file:  f  10Srp0x 

rout  ine  :  al  1  done 

This  will  cause  frames  b-14.  lb-40»  and  44-40  to  be  extracted  from 
the  file  “fl05rp0x“  and  to  be  put  into  the  file  "f0casel*“. 


Example  4: 

It  is  often  best  to  input  everything  on  the  execution  line  tif 
possible! . 

user:  frex  f 105rq0x  1  t  4  box  b44  xl-4lLF!end  /  1.1 

routine  :  Processing  file:  f  105rq0x 

routine  :  all  done 


Example  4: 

This  shows  how  to  obtain  both  frames  with  and  without  alphanumeric 
characters. 

user:  frex  f 103rq0x  box  b44  char+-  /  1.1 
rout ine/user :  /I  t  4  infile  quit  nochar  1  l  4  end 
routine  :  Processing  file:  f 103rq0x 

routine  :  Processing  file:  f 10Srq0x 

routine  :  all  done 
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In  the  output  file  " f0char+- the  first  three  frames  are  identical 
to  those  in  "flBSrqBx"!  but  all  alphanumeric  characters  in  frames  4 
to  b  are  eliminated. 

Warning!  some  tmaybe  all!  characters  generated  by  DISSPLh  are 
draun  and  hence  cannot  be  eliminated. 


Example  5! 


This  is  an  example  of  having  a  file  family  as  input. 


user 

rout  ine/user 
rout ine 
rout ine 
rout ine 
rout ine 
rout ine 


frex  family  fU35arBl  box  b 44  ar/'fam  /  1  .1 

>{  t  5 m  end 

Processing  file!  f  IkJSarBl 

Next  file  in  family!  f  U35arB4 
Last  file  in  family:  f  105ar(3x 
UhRNINL!  INPUT  FILE  CONTAINS  DNLr  444  FAMES 
all  done 


The. warning  given  above  is  to  tell  the  user  the  number  of 
frames  extracted  when  reaching  the  end  of  file. 
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(2)  COMOUT 

Cinco  IliTOUT  creates  a  header  end  a  trailer  for  every  fllo>  it  is 
uosirabl o  to  combine  severe’  fiies  into  one  before  using  NETOUT. 

COMOUT  cn.ibinos  any  number  of  files  into  orio  file;  yot  cacn  file  con 
li.ivo  it:;  own  options.  COMOUT  is  available  on  both  the  CRAY--!  and  the 
C.\C~rijG0.  To  got  it,  typo’, 
i )  Ur.  th..*  Lif'nY—1  — 

rf  ivr.i  road  1222  .Cray  comout (LF) end  /  tv 

>. 1  r 'r.  •  lie  Cl  "i'uGO 

i  i :  i>r.i  rs.  at!  1222  . sivjTG  c  i  rr, out  CLF lend  /tv 

This  document  corresponds  to  the  CQMulJT  version  of  August  10,  1979. 
La  tei  VO!  .ions  of  f.'UliOU  will  bo  stored  in  i-'ILEM  directories  .cruy 
ant!  .sc:>-i. (for  !T:.:Y  .a  VCCD  versions,  respectively)  of  user 
number  ’lin?  rr.  cr  should  periodically  choolt  the  date  o  <  these 

;  (»  •  U-..I  lic-!j  12.'.;  .rrc;y  r  omort  tLlr)!iou  1222  .say7o  comout)  to  see 

Mio  j.r  • . ..'r.!,i:;.  have  beoii  updated.  The  rile  COMOUT  (on  both  machine:,) 
«.i  i.ILi  v  1  lor.  esi.tuinn  the  icte.i  sources  cis  well  as  the  latest 
'..-■i ioi i  o *  this  dorur.urrtcit ioi>.  To  obtain  this  document,  typo; 
lib  ;o;iioult:.i-'>:c  coiwt/dse. (LI*  > end  /  t  v 

-m  FuscLi  omo/duc  bo;  b..r.  comout  /  t  v 


There  are  on  -  i.-..  commands  ii,  COMOUT?  they  are: 

global s:  -  Thu  rest  of  the  input  lino  is  assumed  to  be  NETOUT 
••  ,;•••  ions  fer  the  combineci  file,  fit  most  5  options 
ril  f.  bo  i...-...ccpto.  i>  "global s s "  occurred  more  than 
i.i.'  •  c t ;  or  overrides  the  former. 

end  -  Terminates  COMOUT 


(1)  Doth  cc;, /.rinds  hove  to  be  the  first  symbol  of  ci  input  line. 

(2)  The  user  eon  enter  cs  many  lines  as  ho/she  wants.  Each 
input  lino  that  does  not  contain  a  command  has  to  be 
veil  id  for  NETOUT.  Lino  (cod  anc!  escape  Keys  are  not 
allowed.  The  prompt  from  COMOUT  is  “>*. 

Destination  (silo)  and  bo:;C-.id  should  be  entered  from  input 
lines  oilier  then  that  global  options  are  specified. 
Destination  must  ho  the  first  symbol  of  a  input  line,  and 
Ceroid  must  be  the  lust  three  symbols  of  a  input  line. 

If  they  are  specified  more  than  once,  only  the  f irst  one 
is  t ebon. 

(4)  “global s:  ",  destination,  and  box&id  are  all  optional. 

<S)  Cn  the  .COO,  the  "C  “  option  is  specified  for  "global s: " 
internal ! y.  ftny  option  specified  under  “global s: *  that 
conflicts  the  “0."  option  is  not  allowed. 


Example: 

user 

rout ino/n,or 
ru'd  inn/ 1 .  :er 
i  uu  I  ill.  ■  II.  ..  I 
r  out  it  u* 

t  Ol  d  ii'iu 


comout  global 5 :  uc. 

'■ .'onio, 'doc  box  1.22  comout 
.'I’/h  ('iidwilli.  /:; 

.  ei  !■  I 

'•ii.. nodi'..  :  fill)-  id  is  comout 
nl )  dune 


-  rxucb/Ka 
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(3)  DB 


DS  is  ci  simple  dale  Peso  mcncgsmeri l  system  that  uses  a  hashing 
techanisr.i  to  store  and  retrieve  a  record.  D3  is  available  cn 
he  CSC-FouS  and  cc.n  be  obtained  by  typincj: 
tiler, i  read  1222  .scyFS  c'b(LF)snci  /  tv 


Tnis  report  correspond ing  to  the  DB  version  of  October  2. 

13F9.  Later  version  of  D3  will  be  stored  in  rILZH  directory 
.scy?3  of  user  number  1222.  The  user  should  periodically  ehecK 
the  date  of  this  file  Cfi’.em  !’,os  1222  .scyFS  db)  to  see  if  the 
program  has  boon  updated.  The  file  D3  is  a  L13  file;  it  contains 
the  latest  sources  as  well  as  the  latest  documentation.  Tb 
obtain  the  documental  ion.  type: 
lib  db(LF>:<  db/docCLFJond  /  tv 

nelout  Huso]  dbAloc  uic.  Lswp.0  box  nnn  db  /  t  v 


1-1 .  I  he  ACCESS  option 

The  ACCESS  option  is  used  to  create  and/or  to  update  the  data  !:cse. 
It  has  the  following  syntax: 
db  access  <neu  or  old>=<c:cta  base  file  name>  /  i  v 

If  HEU  is  specified,  the  user  will  be  ashc-d  to  define  the  data 
'rase.  The  following  are  the  prompted  messages  and  their  meanings: 

field  name:  -  Each  record  is  composed  by  fields.  A  maximum 

of  20  fields  is  allowed.  Eccn  field  name  cen 
have  up  to  £3  characters  or  5  words  Cwith  each 
word  of  less  titan  ID  characters) . 

Mote  that  the  first  field  is  always  used 
as  the  hash  Key;  therefore  the  contents 
of  that  field  has  to  be  unique . 
field  length:  -  The  length  in  characters  of  the  field.  The 

maximum  length  allowed  is  £3  shared ors. 
new  field  name:  -  This  will  appear  only  if  the  last  f.oiu  name 

entered  exists. 

This  process  can  be  terminated  by  entering  "end'1  to  ins  prompt 
1  f ield  name : " . 

After  the  data  base  is  defined.,  cr  if  OLD  is  specified,  the  prompt 
wilt  appear  and  the  user  ccn  then  update  the  data  base  by  using 
the  foil  lowing  commands: 

add  <nash  l;oy>  -  Ac'd  a  record  with  the  first  field  =  <hcsh  Kcy>. 

del el e  thash  Kcy>  -  Delete  the  record  with  the  first  field  = 

<h ash  Uoy>. 

update  lhash  l;ey>  -  update  the  record  with  the  first  field  ° 

<has!i  i.sy>.  To  chcnco  ir.e  first  field,  the 
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rhoui  shcsh  hoy> 


user  has  to  delete  the  tchole  record  end  add 
ci  new  one.  The  user  may  nit  on.y  ino  rot  urn 
Key  to  those  fields  that  remain  unchanged. 
Print  the  record  uiih  the  first  fioiu  - 
<hcsh  Key > . 


done 

"srm  inc'-.G 

D3. 

ample: 

user 

db  access  r.ou= 

friends 

rout ins/usor 

f  i  e  1  cl 

ncrr.3 : 

name 

rout ine /user 

f  i  o  l  d 

1  sngth : 

/.3 

rout ine/user 

fie!  c: 

nan*.  3; 

phene  r.u;i-.bor 

rout  ine/user 

field 

length: 

12 

rout ine/user 

f  iold 

r.crna : 

street 

rout ine/user 

field 

length: 

£3 

rout ine/user 

f  iold 

r.crns : 

City 

rout ine/user 

field 

length: 

25 

rout ine/user 

f  1  e  1  a 

name : 

stale 

rout  ine/user 

f  ;  Q 1  d 

1 ongth : 

G 

rout  ir.o/user 

t  ieid 

narns : 

SiO 

rout ine/user 

field 

length: 

5 

rout ine/user 

f  ieid 

nems : 

end 

rout ine/user 

'add  ftu-Yeuncj, 

II.  Stephen 

rout  ir.o/user 

phone 

number : 

41S-C‘-,.;-wC,i- 1 

rout ine/user 

st  roe 

.  ; 

Internet  icncil 

rout ine/user 

city: 

EorKoloi: 

rout ino/usor 

state 

Cf) 

rout ine/user 
rout ine/user 

sip: 

>  clone 

S4F23 

rout ine 

ail 

:!cno 

It  is  i [-.-.portent  to  note  that  c:n  upper  case  chcract 
internally  toy  three  characters  -  CO  plus  the  Icv.-er 


:resencs 
;r  acner. 


3.  The  PRINT  cp-ticn 

To  generate  a  list  of  the  data  Sets®,  cr.e  uses  the  PRINT  option 

a s  ( c ! i quo ■ 

ab  print  sclata  bc.se  file  ncme>  [output  tile  name  [si  tool  /  t  v 

if  the  output  file  name  is  omitted,  it  defaults  to  "dbcut.p". 
if  site  is  given  after  the  output  file  name.  the  output  uill  be 
sent  to  the  destination  <sito>  by  running  the  NITC'JY  pr ogrc;..i. 

The  user  uill  then  be  ashed 'the  fields  that  he  uculd  1 iho  Is  be 
sorted.  This  should  be  entered  from  the  highest  priority  is  the 
louest  and  terminatc-d  by  entering  "end". 


3.  'ample: 

user 

rout ine/user 
rout ine/user 
rou  t  ine/user 
rout  ine 
rout ine 


a'b  print  friend  friend. p 
field  to  be  sorted:  state 


field  to  be  sorted t  name 
field  to  be  sorted:  end 


cor  i  ;cr!n;  r.c;  i  cci 


rr.zil .  u 


all  der.o 
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(k)  TIME 


TIME  allows  any  program  (except  TIME  itself)  to  be  run  under 
its  control.  It  prints  the  time  used  by  the  controlleo  when  the 
controllee  terminates.  TIl'iE  is  available  on  the  C3rtY-l  and  can 
be  obtained  by  typing! 

rfilem  read  1222  .Cray  time(LF)ond  /  i  v 


This  document  corresponds  to  the  TIME  version  of  July  23.  197S. 
Later  version  of  TIME  uiili  be  stored  in  rILEM  directory  . Cray  of 
user  number  1222.  The  user  should  periodically  checK  the  date  of 
this  file  (filern  how  1222  .cray  time)  to  see  if  the  program  has 
been  updated.  The  file  TIME  is  a  LIB  file;  it  contains  the  latest 
sources  as  well  as  the  latest  documentation.  To  obtain  the 
docurnentcit  icn.  type: 

i  ib  t  ime CLF) x  t  ime/dec(LF)ond  /tv 
net out  Cusc]  lime/doc  box  bnn  time  /  tv 


There  is  only  cno  command  in  TIME  -  fin.  to  terminate  the  procjrcm. 
The  user  can  sequent icl 1 y  run  as  many  programs  (control  1 ees)  as 
he/ she  wishes.  Before  entering  the  name  of  the  next  controllee. 
the  user  should  wait  for  the  prompt  “»"  to  appear. 
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K.  REPORT  GENERATION 

Alex  Friedman 

Recent  developments  have  made  report  generation  on  the  M/FECC 
computer  facilities  more  convenient  for  users  at  sites  serviced  by 
mini-USC's,  such  as  U.C.  Berkeley.  The  first  of  these  Is  the 
availability  of  the  output  options  "nice"  and  "xnice"  in  the  "PRINT" 

t 

command  of  TRIX  AC;  these  options  send  output  to  the  NIPS  printer  with 
drawn  type  fonts.  To  find  out  more  about  them,  obtain  a  copy  of  the 
latest  RED  report,  as  described  in  the  previous  OPR.  Note  that  MIPS 
output  arrives  more  promptly  than  high-quality  hardcopy  output,  and  is 
far  less  costly.  Also  note  that  plot  files  generated  in  any  manner  can 
be  sent  to  the  NIPS  by  use  of  utility  routine  "NIPIT”,  and  so  if  the 
user  prefers  to  invoke  REDPP  directly  rother  than  through  the  TRIX  PRINT 
command,  this  is  possible. 

The  second  development  is  the  imp! ementat ion  of  the  "TURN."  option 
of  NETPlOT,  the  routine  which  is  used  to  send  graphics  files  to  the 
local  Versatec  pr inter /pi ot ter .  This  option  is  useful  for  all  plotted 
output,  as  the  plots  are  rotated  by  90  degrees  and  so  the  user  can  leaf 
through  the  output  in  book  form.  However,  it  is  particularly  useful  for 
plotting  text,  as  one  can  now  obtain  output  where  the  pages  cppear  in 
"forward"  order.  The  "ampersand  zq"  command  to  REDPP  can  be  deleted.  A 
sample  COSMOS  deck  to  perform  the  appropriate  processing  on  the  file 
"example"  of  the  previous  OPR  is  the  revised  file  "swapexcm",  which 
contains: 
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•select  p=swoprec 

♦netout  ucb  example  s,  ulc.  box  b34  exomple 

»tri*  oc 

•-red? 

*  -a 

••  -examp  I  e 

«-r I  1,2 

•-/popersize  66  72  .04  .1  .32  0.  0. 

«-nf ( tempO) 

« - f  o  r ma  t ( t  empO . t  emp 1 ) 

»  -end 

•redpp  tempi  startp.  1  pcperr.  keep,  defont.  2  2hi  ts.  & 
•-filesize.  1360000  frSCbrk.  1300000 
»  — box  b34  example 

♦netplot  ucb  alwith.  fOhcy  fxhcy  f.  I.  turn,  box  b34  example 
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Station  Ill 

PLASMA  SIMULATION  TEXT 

One  publisher  has  manuscript  for  review,  looking  favorable.  As  this 
has  dra^ea  on  (their  problem),  up-dating  will  be  required,  when  and  if  accepted. 
We  have  been  out  of  copies  for  months  and  do  not  plan  to  reprint,  at  present. 

Station  IV 

SUMMARY  of  REPORTS,  TALKS,  PUBLICATIONS  IN  PAST  QUARTER 

Abstracts  follow  of  papers  submitted  (this  quarter)  for  Division  of 
Plasma  Physics  Meeting,  American  Physical  Society,  Boston,  Mass.,  Nov.  12-16, 

1979  (next  quarter) . 

Jae  Koo  Lee  and  C.  K.  Birdsall,  "Velocity  Space  Ring-Plasma  Instab¬ 
ility,  Magnetized,  Part  I:  Theory,  Part  II:  Simulation",  Phys.  FI.  22, 

7,  pp.  1306-13^,  1315-1322,  July,  1979. 

H.  Stephen  Au-Yeung  and  Alex  Friedman,  "Solver:  An  Analytic  Function 
Root  Solving  and  Plotting  Package",  ERL  Memo  No.  UCB/ERL  M79/55,  31  August 
1979. 


Saturation  of  the  Lower-Hybrid  Drift  Instability. 
YU-JIUAN  CHEN  and  C.K.  BIRDSALL,  U.C.  Berkeley*— The 
linear  properties  and  the  saturation  mechanism  of  the 
lower-hybrid  drift  instability  are  studied  using  a  ID 
particle-hybrid  simulation.  The  model  is  a  slab  with  a 
constant  density  gradient;  the  ions  are  unmagnetized  par¬ 
ticles,  shielded  by  the  strongly  magnetized  electrons 
through  the  linear  electron  susceptibility,  Xe-  Ions  are 
initially  in  a  steady  equilibrium  state  with  the  ion  dia¬ 
magnetic  drift  velocity  cancelled  by  the  E*B  drift,  cor¬ 
responding  to  electrostatically  confined  ions.  At  small 
amplitudes,  the  simulation  shows  good  agreement  with  lin¬ 
ear  theory,  such  as  the  linear  growth  rate,  the  real  fre¬ 
quency,  and  the  influence  of  finite  beta  effects  associ¬ 
ated  with  the  nonresonant  ?B0  electron  orbit  modifications. 
At  large  amplitude,  allowing  only  a  single  mode,  it  is 
found  that  the  end  of  wave  growth  is  due  to  ion  trapping, 
even  when  a  wide  band  (Au/vy)  of  the  mode  occurs,  with  the 
growth  rate  y  comparable  to  the  wave  frequency.  Contrast 
with  the  end  of  growth  by  quasilinear  diffusion  will  be 
given. 

Particle  Simulation  of  Instabilities  due  to  Steep 
Density  Gradients.  JAE  K00  LEEtandC.K.  BIRDSALL,  U.C. 
Berke  1  ey* — Instabilities  may  occur  in  co  1  1  i  s  i  on  less  Max¬ 
wellian  magnetized  plasmas,  driven  by  the  free  energy  as¬ 
sociated  with  a  spatial  density  gradient.  Electrostatic 
particle  simulations  were  used  to  study  such  instabilities 
with  both  species  magnetized  and  treated  fully  nonlinearly. 
During  the  linear  stage,  simulations  showed  exponential 
growth  in  time  with  the  growth  rates  in  fair  agreement 
with  a  linear  nonlocal  theory,  while  the  real  parts  of 
frequencies  were  not  well  resolved  in  the  short  growth 
period.  The  simulation  saturation  levels  were  somewhat 
above  those  predicted  by  existing  nonlinear  theories.  At 
the  time  of  saturation,  the  phase  space  pictures  of  elec¬ 
trons  and  ions  show  bunching  i n  some  cases.  Clear  dis¬ 
tinction  between  drift  cyclotron  and  lower  hybrid  drift 
instabilities  was  not  possible;  both  may  have  been  pre¬ 
sen  t . 

Field  Reversed  Ion  Ring  Stability  -  Recent  Results; 
Erqodic  Orbits  and  Particle  Simulat ion.*  A.  FRIEDMAN. 

U.  C.  Berkeley.  J.  DENAVIT.  Northwestern  Univ..  and  R.  N. 
SUDAN.  Cornell  Univ.  Ue  present  new  results  regarding 
stability  of  a  f iel d-reversed  ion  ring  in  a  dense  plasma, 
obtained  by  numerical  simulation  using  RINGHYBRID.  a 
linearized  3D  hybrid  code1.  The  ring  is  moderately  thick, 
with  effective  aspect  ratio  of  order  4:1.  and  reversal 
factor  1.35  on  axis.  Nonax isymmetr ic  modes  of  azimuthal 
number  ?  are  studied.  The£=l  redial  mode  (precession)  is 
stable,  whereas  the  axial  mode  (tilt)  is  unstable.  Axial 
kink  modes  with  St  >2  and  radial  kink  modes  with  X >3  are 
stable,  as  predicted  on  the  basis  of  thin-ring  theory1'3  . 

Ue  have  observed  effects  of  ergodic  particle  orbits 
in  our  simulations.  In  the  nonl inear  2D3V  zero  order  be¬ 
havior  the  main  effect  is  a  loss  of  left-right  symmetry 
due  to  the  exponential  divergence  of  “neighboring*  mirror 
image  tra jector ies .  However,  in  linearized  codes  where  the 
perturbat ion  quantities  represent  displacements  between 
neighboring  orbits,  the  collective  behavior  can  be  masked 
in  some  cases  by  rapid  single-particle  orbit  separation. 

*5upp.by  USDOE  Conlrs.  DE-AS03-763F0G034xDE-AT33-76ET53064 
and  EY7S-S-02.2208. 

'A. Friedman.  R.N. Sudan.  J.Denavit.  Cornell  Lab.  of  Plasma 
Studies  F.ept.#263  (1979).  submitted  to  J.  Comput .  Phys. 
1R.V. Love lace.  Phys.  Fluids  19.  723  (1976). 
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